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Abstract 

We review recent developments in the modelling of the 
phase diagram and the kinetics of crystaUization of car- 
bon. In particular, we show that a particular class of 
bond-order potentials (the so-called LCBOP models) 
account well for many of the known structural and ther- 
modynamic properties of carbon at high pressures and 
temperatures. We discuss the LCBOP models in some 
detail. In addition, we briefly review the "history" of 
experimental and theoretical studies of the phase be- 
haviour of carbon. Using a well-tested version of the 
LCBOP model (viz. LCBOPI+) we address some of 
the more controversial hypotheses concerning the phase 
behaviour of carbon, in particular: the suggestion that 
liquid carbon can exist in two phases separated by a 
first-order phase transition and the conjecture that di- 
amonds could have formed by homogeneous nucleation 
in Uranus and Neptune. 

1 Introduction 

Carbon exhibits a rich variety of solid structures: Some 
are thermodynamically stable, most are not. To be 
specific: solid carbon can be found in the two well- 
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known crystalline phases, diamond and graphite, and 
in amorphous states, such as glassy carbon and car- 
bon black. Furthermore, the existence of additional 
metastable solid phases at relatively low pressure, the 
so-called carbynes, is still hotly debated [Ull]). In ad- 
dition to the bulk phases, there are the more recently 
discovered fuUerenes, Cgo and C70 [3], nanotubes [4], 
and graphene [5]. 

The reason why a simple element such as carbon can 
manifest itself in so many different forms is related to 
its unusual chemical properties: carbon exhibits three 
different possibilities for covalent bond formation: sp^ 
hybridization appears in diamond, sp2 hybridization is 
found in graphite, graphene, nanotubes, and fullerenes, 
whilst in carbynes, C should exhibit sp hybridization. 
Because of their high cohesive energies and concomi- 
tant high activation energies that must be overcome in 
structural phase transformations, carbon polymorphs 
often exist in metastable form well inside pressure- 
temperature regions where another solid form is ther- 
modynamically stable. For example, it is well known 
that diamonds survive at normal P-T conditions, where 
graphite is the thermodynamically stable phase. Con- 
versely, graphite tends to persist at very high pressures, 
deep into the diamond stability region of the phase di- 
agram. 

It is also interesting that, at zero pressure and tem- 
perature, graphite and diamond have a very similar 
(and quite large) binding energy per atom, i.e. 7.37 
eV (graphite) vs 7.35 eV (diamond). This fact might 
suggest (and it has indeed been suggested) that also in 



disordered phases like the hquid, the two local struc- 
tures - graphite- like and diamond-like ~ could compete. 
In fact, as we discuss below, the possibility of the ex- 
istence of two distinct and partially immiscible liquid 
phases of carbon has been a subject of much debate. 

In a liquid-liquid phase transition (LLPT), a liquid 
substance displays an abrupt change in some local or 
global property within a narrow band of pressures and 
temperatures. Local properties that may change in 
a LLPT are the local coordination or hybridization, 
typical global properties that are affected are the den- 
sity or the resistivity. LLPT's in dense, atomic liq- 
uids are typically difficult to probe experimentally: the 
candidate transitions often occur at extreme pressures 
and/or temperatures or appear in metastable regions of 
the phase diagram (and may be hidden by competing 
solidification). Evidence for LLPT's have been found 
for a number of atomic systems, such as Cs [6], As [7], 
Bi [5], Ge 0, Hg [in], S [H], Sb [H], Se [H], Si [Tmi], 
Sn [IS], H2 [n], I2 [H], N2[Ill2n]. The best estab- 
lished experimental example of a LLPT in an atomic 
liquid is the case of phosphorus. A transition between 
a fluid of tetrahedral P4 molecules and a network form- 
ing (and metallic) liquid was predicted on theoretical 
grounds [HJ [25] , and subsequently verified experimen- 
tally [331 [21] • The LLPT in phosphorus has been anal- 
ysed in several numerical studies [3S1 HH] HTJ JH] ■ Many 
other network-forming liquids are also expected to ex- 
hibit LLPT'a: first and foremost water [29l [30], but 
also Si02 [31] and Ge02 [32]. Although considerable 
progress has been made in the theoretical description 
of LLPT's [331 [3113313a 133 131] , a unified theoretical 
picture is still lacking. 

In this review we discuss the phase diagram of solid 
and liquid carbon at high pressures and temperatures 
on the basis of the results of numerical simulations; 
both quantum and classical. We present evidence that 
the presence of graphite-like and diamond-like local 
structures in the liquid does not give rise to liquid- 
liquid demixing but that the predominant local struc- 
ture in the liquid varies strongly with pressure. The 
fact that the liquid is locally either graphite-like or 
diamond-like has dramatic consequences for the nucle- 
ation of the diamond phase, a finding that may have 
some consequences for our understanding of carbon- 
rich planets or stars. Wherever possible, we discuss 
our own results in the context of the relevant litera- 
ture about the carbon phase diagram, about a possible 
LLPT in this system and about the possibility of dia- 
mond formation in planetary interiors. 

In section [2] we give an overview of the bond-order 
potential ("LCBOP") that was used to compute both 
the equilibrium phase diagram of carbon and the path- 



way for diamond nucleation in liquid carbon. We also 
discuss in some detail the different variants of the 
LCBOP potential [3S1 HO] and explain the rationale 
behind the choice of the present LCBOP potential. 

In section [3] we briefly summarize some of the earlier 
ideas about the phase diagram of carbon (in particular, 
about the crystalline phases and the liquid). We pay 
special attention to the slope of the diamond melting 
curve and the (possible) heating-rate dependence of the 
graphite melting curve. 

In section 13.21 we report our results concerning the 
phase diagram, in the context of recent first-principle 
simulations and experiments. 

In section 14.11 we review the arguments that have 
been put forward to support the idea that liquid carbon 
can undergo a LLPT. We argue that, to the extent 
that we can trust the present models of liquid carbon, 
a LLPT in carbon is not be expected. 

In section[5]we discuss our numerical results concern- 
ing the (homogeneous) nucleation of diamond from the 
bulk liquid. In particular, we discuss in some detail the 
numerical approach that was used in Ref. 41]. In ad- 
dition, we focus on the structural analysis of the small 
solid clusters in the liquid and we discuss the impli- 
cations of our findings for the formation of diamonds 
in carbon-rich star systems and the interior of giant 
planets. 

2 The LCBOP-family 

Whilst the crystalline phases of carbon can be simu- 
lated by traditional force fields that do not allow co- 
ordination changes, a study of the liquid phase and, a 
fortiori, of phase transformations between phases with 
different local coordination, requires a potential that 
can describe carbon in different coordination states. 
The LCBOP potential was designed with this objec- 
tive in mind. LCBOP stands for "Long range Car- 
bon Bond Order Potential", and represents a bond- 
order potential for pure carbon that includes long- 
range (LR) dispersive and repulsive interactions '39' 
from the outset. We stress that LCBOP is not based 
on an existing short-range (SR) bond-order potential 
to which LR interactions have been added a posteriori, 
although such an approach has been proposed in the 
literature [121 US] HI] . The latter approach requires a 
rather special procedure to avoid interference with the 
SR potential and suffers from a loss of accuracy. In the 
LCBOP these problems are circumvented in a natural 
way. 

We note that the term "long range" may be confus- 
ing in this case, as the cut-off of the LR potential in 
LCBOP is only 6 A. Usually, the term "long range" is 



only used for interactions with a much longer range, 
such as Coulomb interactions. Here we use it as a 
synonym for "non-bonded" , referring to a range much 
larger than the typical distances between chemically 
bonded atoms. 

After the introduction of LCBOP in Ref. [3^], a 
number of significant modifications have been intro- 
duced in order to improve its description of all carbon 
phases, including liquid carbon. To facilitate the dis- 
tinction between the different LCBOP potentials, the 
various versions have been named LCBOPI, LCBOPI+ 
and LCBOPII. LCBOPI, introduced as LCBOP in Ref. 
[39], does not include torsion interactions. As tor- 
sion interactions were shown to play an important role 
in liquid carbon [¥5], we introduced a refinement of 
LCBOPI, called LCBOPI"*", when we performed our 
first study of liquid carbon [IS]. LCBOPI"*", includes, 
among other changes, conjugation dependent torsional 
interactions. Clearly, describing the liquid phase re- 
quires a robust form of the potential in order to deal 
with configurations that are quite unlike the regular 
topologies in crystal lattices. LCBOPII addresses this 
problem: it includes several important improvements 
over LCBOPI"*". An important innovation in LCBOPII 
is the addition of so-called middle range (MR) interac- 
tions, introduced to bridge the gap between the extent 
of the tail of the covalent interactions as found in ab- 
initio calculations, (up to 4.5 A in certain cases) and 
the rather short cut-off of only 2.2 A in LCBOPI+ for 
these interactions. 

In the remainder of this section we give a brief step- 
by-step description of the transition from bond-order 
potentials (BOPs) to LCBOPI+. All the results pre- 
sented in later sections, concerning the phase diagram, 
the liquid structure, and the nucleation issues are ob- 
tained with this version of the potential. In Appendix 
A we discuss the LCBOPII potential, with a short ac- 
count of results obtained with this refined version of 
the potential. We aim at giving the flavour of the po- 
tentials in a mainly descriptive way with graphical il- 
lustrations, minimizing mathematical formulation. 

2.1 Bond-order potentials 

A bond-order potential (BOP) is a reactive potential, 
i.e. able to deal with variable coordination. It provides 
a quantitative description of the simple idea that the 
bonds of an atom with many neighbours are weaker 
than those of an atom with few neighbours, as the co- 
hesive ability of the available electrons has to be shared 
among the neighbours. For carbon, each atom delivers 
four valence electrons. If these four electrons have to 
make the six bonds in a simple cubic lattice, then it is 



evident that each of these bonds is weaker than a bond 
in the diamond or graphite lattice with coordinations 
4 and 3 respectively. On the other hand, the number 
of bonds is larger for the simple cubic lattice. So there 
is a balance to be made, which in the case of carbon 
has the result that graphite is the most stable phase at 
ambient pressure. 

A thorough analysis of these bonding properties, 
based on a quantum mechanical description, has been 
given by Anderson [47l|48l|49] and AbcU [50]. In their 
description, that forms the basis of the tight binding 
models, the electronic wave function is approximated 
as a sum of localized atomic orbitals. Abell showed 
that for a regular lattice, i.e. with an identical environ- 
ment for each atom, and within the assumption that 
the overlap integral for orbitals on different atoms is 
non- vanishing only for nearest neighbours, the binding 
energy per atom is given by: 



E, = ^ZiqVn{r)+bVAir)) 



(1) 



where Z is the number of nearest neighbours, q is the 
number of valence electrons per atom, V^{r) is a two- 
body potential describing the core repulsion, VA{r) is 
a two-body attractive potential, and b is the so-called 
bond order, a many-body term dependent on the lo- 
cal environment of the atoms. Abell also showed that 
the coordination dependence of b is fairly well approx- 
imated by: 



b = b{q,Z) = a{q)Z-^^^ 



(2) 



with a{q) a function of q, specified in Ref. 3D]. As- 
suming exponential functions VR(r) = Aexp{—6r) and 
Va(^) — — Sexp(— Ar), with A, B, A, and 9 fitting 
parameters, as a reasonable approximation for over- 
lapping atomic orbitals from atoms at distance r, and 
defining S — 9/X, some algebra leads to a total binding 
energy given by: 
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and an equilibrium nearest neighbour distance given 
by: 
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where C and C" are constants. Eq. |3| implies that for 
S > 2 high coordination structures (close packing) are 
favoured (metals), whereas for S < 2 the dimer will 



be the most stable structure (extreme case of covalent 
bonding). As, in general, the repulsion falls off (much) 
faster than the attraction, i.e. 9 > X {S > 1), Eq. [3] 
implies that Veq is monotonically increasing with coor- 
dination. Combining Eqs. [3] and |4] yields a simple rela- 



tion between Eb and req, namely Eb esc exp {9 — 2A)r, 



eq- 



value 4 of carbon. The character of a certain bond ij, 
and its conjugation number A'^"'"-', a number between 
and 1 quantifying effects beyond nearest neighbours, 
is determined by the sum of the electrons supplied by 
atom i and atom j. The interpolation model is further 
illustrated in Fig. [H 



In principle BOPs are based on the above bonding 
ideas. However, the transferability to different types 
of structures and materials has been greatly enhanced 
by a quite reasonable extension in the functional form 

of the bond order b. The simplest bond order bij for a 03 LCBOPI"^ 
bond ij according to a BOP in the style of Tersoff [5T| 
and Brenner [52] reads: 



bt] = 
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where the sum runs over the nearest neighbours other 
than j of atom i, G{9ijk) is an adjustable function of 
the bond angles 9ijk and e is a negative exponent but 
not necessarily -1/2. Taking a constant G{9ijk) — 1 
and e = —1/2 yields b ~ aZ^^'^, i.e. one recovers 
the functional form Abell found for b (Eq. [5|). This 
form includes the effect of bond angles in a natural 
way, and has the ability to fit a large set of data quite 
well, explaining the success of BOPs. 



The potential LCBOPI+ is given by LCBOPI supple- 
mented with torsion interactions and a correction of 
the angle dependent part of the bond order for con- 
figurations involving low coordinations and small an- 
gles. Similar modifications were also included in the 
REBO potential ^K3], although the torsion term con- 
tains a significant difference. For LCBOPI+, following 
the results of ab-initio calculations, the shape of the 
torsion energy curve as a function of the torsion angle 
depends on conjugation, whereas for the REBO poten- 
tial the curves for a double and a graphitic bond are 
equally shaped but scaled (see top panel of Fig. [20|) . 
Details on LCBOPI+ are given in the appendix A of 
Ref. [H] andinRef. [55]. 



2.2 LCBOPI 

The main innovative feature of LCBOPI [39| concerns 
the treatment of the LR van der Waals interactions. 
One of the challenges here is to add LR interactions 
which describe not only the interlayer graphitic bind- 
ing but also the rather strong 7r-bond repulsion in 
graphite for decreasing interlayer distance without pay- 
ing a price in the accuracy of the covalent binding 
properties. A correct description of these interactions 
requires a LR potential, V^'''(r) that is repulsive in 
the distance range corresponding to the second nearest 
neighbours (in diamond and graphite). In the LCBOP- 
family, to get the right equilibrium lattice parameter, 
this extra repulsion has been compensated by a some- 
what stronger attractive part of the covalent interac- 
tion, achieved by an appropriate parametrisation of the 
SR potential. This is schematically illustrated in Fig. 

HI 

Another feature of LCBOPI is that it contains a rea- 
sonable, physically motivated interpolation scheme for 
the conjugation term to account for a mixed saturated 
and unsaturated environment. In this approach each 
atom supplies a number of electrons to each of its bonds 
with neighbouring atoms according to a certain distri- 
bution rule, the total sum being equal to the valence 




Figure 1: Schematic representation of the LCBOP ap- 
proach to the inclusion of LR interactions. To preserve 
the right equilibrium bond distance for a given bond 
ij, the repulsion of atom i due to the added LR inter- 
actions with the atoms h and I2 (represented by the 
force fii = fjij^ -I- fi/j) is compensated by a stronger at- 
tractive part in the SR interaction between i and j, 
represented by the extra force Afij. The same holds 
for atom j. For convenience, here only the main, re- 
pulsive LR interactions are drawn, whereas in reality 
this compensation is for the sum of all LR interactions. 




Figure 2: Schematic presentation of the interpolation 
approach in the determination of A^™"-' for mixed co- 
ordination situations (c and d). Black atoms are satu- 
rated atoms, white atoms are unsaturated. The num- 
bers near the atoms i and j indicate the distribution 
of electrons among their bonds. Each single bond, i.e. 
a bond with a saturated neighbour, takes 1 electron. 
The remainder of the electrons (in total 4) are equally 
shared among the bonds with unsaturated atoms. A 
linear dependence of N^""'-' on the total number of elec- 
trons involved in the bond ij with < Nf°^-' < 1 as 
an additional constraint leads to the given values for 

fqconj 



The phase diagram of carbon 
at very high pressures and 
temperatures 



In this section we give a review of experimental and 
theoretical works aimed at determining the phase be- 
haviour of carbon at high temperatures and pressure. 
We follow a "historical" approach, starting from the 
beginning of the twentieth century, up to the most re- 
cent results coming from experiments and computer 
simulations. A historical approach may give a bet- 
ter understanding why certain issues have been, and 
in some case still are, controversial. After setting the 
stage (with a particular attention to the ideas about 
the sign of the slope of the diamond melting curve and 
the long debated issue of the position and nature of the 
graphite melting curve), we focus on the topic of the 
LLPT for carbon. 



3.1 The history of carbon phase dia- 
gram 

One of the earliest phase diagrams of carbon appeared 
at the beginning of the twentieth century, and is due 
to H. Bakhuis Roozeboom |56J, who estimated the 
phase behavior of carbon on the basis of thermody- 
namic arguments. Of the two solid phases, diamond 
was recognized to have a slightly greater vapor pres- 
sure at a given temperature. The temperature of the 
graphite/liquid/vapor triple was believed to be around 
3000 K. In 1909 Tamman j57J postulated the existence 
of a region where graphite and diamond are in pseudo- 
equilibrium. The existence of this pseudo-equilibrium 
region was at the basis of the method of synthesizing 
diamond starting from carbon saturated solutions of 
molten iron, silver, or silicates. In 1938, Rossini and 
Jessup j58j of the U.S. Bureau of Standards used ac- 
curate thermodynamic data to estimate that at K 
the lowest pressure at which diamond would be sta- 
ble against graphite is around 1.3 GPa, and around 
2 GPa at 500 K. In 1939, the Russian scientist Leipun- 
skii 59J published a review of the problem of diamond 
synthesis. On the basis of thermodynamic data, he 
suggested that the melting curve of graphite might be 
at about 4000 K, with possibly some increase with 
pressure. This value for the melting of graphite was 
rather well verified the same year by Basset [BO , who 
established the graphite/liquid/ vapor triple point to be 
at about 11 MPa and 4000 K. In that same publica- 
tion. Basset reported on a rather pressure independent 
melting temperature of graphite at ^ 4000 K, from 
atmospheric pressure up to 0.1 GPa. In 1947 Bridg- 
man [ST] addressed the problem of extrapolating the 
graphite/diamond coexistence curve beyond the region 
where it can be estimated from known physical prop- 
erties (4 GPa/ 1200 K). He concluded that there was 
a possibility that at higher temperatures the rate of 
increase of P with T along the curve would decrease. 
This hypothesis was later supported by Liljeblad [62] 
in 1955, while Berman and Simon [63j in the same 
year came to the conclusion that the best extrapolation 
would be a straight line. Experiments that could decide 
this issue were started by Bundy and coworkers in 1954, 
when they accomplished diamond synthesis by activat- 
ing the graphite-to-diamond reaction with the use of 
different solvent-catalyst metals. The relevant exper- 
imental data were published only much later [641 I65j , 
and are compatible with the Bcrman-Simon straight 
line extrapolation. 

Bundy and his group made also extensive experiments 
on graphite melting at pressures much higher than the 
graphite/liquid/vapor triple point. The determination 



of the graphite melting curve is an experimental chal- 
lenge for several reasons. First of all, to reach pressures 
as high as 10 GPa, the sample must be in direct con- 
tact with a solid container and, because the melting 
temperature are so high, this container must be made 
of a material that is as refractory and inert as possible 
(Bundy chose boron nitride, pyrophyllite, MgO and di- 
amond powder). In addition, both the heating of the 
sample and the observations of the high-pressure/high- 
temperature phase must be carried out very rapidly, 
before the wall material can melt or react with the 
carbon sample. The experiments were performed by 
discharging an electrical capacitor through the sample 
(this procedure is known as flash heating) , and by mon- 
itoring the current through, and the voltage across it 
by means of a two-beam oscilloscope. The discharge 
circuit was designed to have energy insertion in the 
sample within a few milliseconds. The interpretation 
of such experiments, is rather sensitive to the assumed 
pressure and temperature dependence of the material 
under study and on the assumption that the pressure 
of the graphite specimen during rapid heating is the 
same as in a quasi-static process. With these assump- 
tions, Bundy' s experiments gave a graphite melting 
curve as shown in Fig. [31 A maximum melting tem- 
perature of about 4600 K was detected in the region of 
6 GPa to 7 GPa. The presence of a region with a neg- 
ative dT/dP along the melting curve indicates that, at 
those pressures, the density of the liquid at the melting 
temperature is greater than that of the solid. 
Interrupting for a moment the historical order of 
events, we note that, throughout the past century, dif- 
ferent experiments located the graphite melting curve 
at rather different temperatures [H^l [SOI IS3 UHl [Ml ISS 
[701 ini [Za [12 [HmS]. At low pressure, the melting 
temperature (T^) was found at values ranging from 
- 3800 to - 5000 K. Asinovskii et al. [TO] pointed out 
the non negligible dependence of graphite T„i on the 
heating rate of the sample. Specifically, heating times 
of the order of 10~^ s [701 IZ3 yielded estimates of 
Tm ~ 4800 - 5000 K; heating times of the order of 
10-3 s ini [74] suggested T„, ~ 4500 - 4600 K; finally 
experiments with heating times of the order of one sec- 
ond jBO] [00] were consistent with the assumption that 
Tra ~ 3800-4000 K. In ref ?6i|, after a thorough discus- 
sion of the experimental methods, the authors recom- 
mended that only data coming from experiments with 
heating time of the order of seconds or more should 
be accepted. This implied that most of the available 
data on graphite melting had to be reconsidered and 
that the question of the position and the shape of the 
melting curve is still open. On the basis of a series of 
laser induced slow heating experiments [76] (i.e. heat- 



ing times of the order of one second), Asinovskii et 
al. proposed the triple point vapor/liquid/graphite at 
~ 4000 K and 0.1 MPa (i.e. atmospheric pressure), 
in clear contradiction to the commonly accepted val- 
ues [73] of - 5000 K and 10 MPa. The next year 
the same authors [77] published results concerning the 
position of the graphite melting curve. With ohmic 
heating of graphite samples at heating rates of about 
100 K/minute, they found T,„ = 3700 K at 0.25 MPa 
(typically, samples melted after one hour of steady 
heating). 

Coming back to Bundy's work, during the experiments 
on graphite melting Bundy and his group also inves- 
tigated the graphitization of diamond by flash-heating 
under pressure. Small diamond crystals where embed- 
ded in the graphite sample, pressurized and then flash- 
heated. Experiments indicated that there is a sharp 
temperature threshold at which the diamond crystals 
completely graphitized. This threshold is a few hun- 
dreds degrees lower than the graphite melting curve. 
Attempts to obtain direct (i.e. without resorting to a 
catalyst material) conversion of graphite into diamond 
by the application of high pressure date back to the 
beginning of the twentieth century. Success came only 
in 1961, when De Carli and Jamieson [7S^ reported the 
formation and retrieval of very small black diamonds 
when samples of low-density polycrystalline graphite 
were shock compressed to pressures of about 30 GPa. 
Later in 1961 Alder and Christian 79J reported re- 
sults on the shock compression of graphite that were 
in substantial agreement with those of De Carli and 
Jamieson. 

Bundy [80] achieved direct conversion of graphite 
into diamond by flash-heating graphite sample in 
a static pressure apparatus, at pressures above the 
graphite/diamond/liquid triple point. The threshold 
temperature of the transformation was found several 
hundred degrees below the melting temperature of the 
graphite, and decreasing at higher pressures. The 
phase transition was revealed by a sharp drop in the 
electrical conductivity of the samples (that were re- 
trieved as pieces of finely polycrystalline black dia- 
mond). 

By linking his own results with earlier experimental 
and theoretical findings, Bundy 80J proposed in Febru- 
ary 1963 a phase diagram of carbon at high pressures 
that is illustrated in Fig. [3l The diamond melting curve 
was believed to have negative slope by analogy with the 
other Group IV elements [52 [52]), and on the basis of 
evidence collected during the experiments of Alder and 
Christian [70]. 

In 1973 Van Vechten f83^ predicted the phase diagram 
of carbon by rescaling the behavior of other Group IV 
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Figure 3: The phase diagram of carbon at high pres- 
sures proposed by Bundy in 1963 [501 ■ 



elements that are experimentally more accessible, us- 
ing the electronegativity as a scale parameter. In 1979 
Grover [M] calculated a phase diagram by using a 
phenomenological equation of state for the description 
of various solid and liquid phases of carbon. He used 
physically motivated approximations for the free ener- 
gies of the various phases, with parameters adjusted to 
match the available data on the equations of state. He 
concluded that, at all pressures, diamond transforms, 
before melting, into a solid metallic phase. 
On the basis of experimental evidence [SS], in 1978 
Whittaker [IJ proposed the existence of a novel crys- 
talline phase for elemental carbon, called "carbyne" . 
The stability region of carbyne is sketched in Fig. 2]) 
and the structure of this phase (though expected in 
different allotropes) is generally that of a chains of al- 
ternated single and triple bonds, i.e. (— C=C— )„, ar- 
ranged in a hexagonal array bundled by dispersion in- 
teractions. The existence of a carbyne form was later 
questioned by Smith and Buseck [2 , who claimed that 
all the experimental evidence could also be accounted 
for by the presence of sheet silicates. This dispute con- 
tinued to this day: the experimental evidence for the 
existence of carbyne is still being debated. 
In recent years, the experimental effort has focused on 
the collection of reliable data at even higher pressures, 
and on the investigation of the properties of the differ- 
ent phases of carbon at high temperatures and pres- 
sures. This challenging task has been faced both with 
experiments and theory. On the experimental side, the 
development of the diamond- anvil cell [6 for high pres- 
sure physics has made it crucial to know the range of 
stability of diamond under extreme conditions. The 
availability of high-energy pulsed laser sources led to 
new tools for heating up samples at very high temper- 
atures (above the graphite melting curve) [SS] . These 



techniques were immediately applied to the determina- 
tion of the properties of liquid carbon (i.e. whether it is 
a conducting metallic liquid or an insulator). Unfortu- 
nately, due to the difficulties in interpreting the results 
of these experiments, the nature of the liquid state of 
carbon is still not characterized experimentally. 
On the theoretical side, the appearance of ever more 
powerful computers made it possible to use electronic 
density-functional (DF) theory [STl [HH] to predict the 
properties of materials under extreme conditions. In 
1983 Yin and Cohen [351 studied the total energy ver- 
sus volume and the free energies versus pressure for 
the six possible lattices of carbon (fee, bcc, hep, sim- 
ple cubic, /3-tin, diamond). The study was carried out 
by using ab initio pseudopotential theory (this permits 
the investigation of the properties of the atomic sys- 
tem at K). Yin and Cohen found that the calcu- 
lated zero-pressure volume for diamond is cither close 
to or even smaller than those of the other five phases. 
This is different from what is observed for the other 
group IV elements. Si and Ge, and defies the common 
notion that diamond is an open structure and should 
have higher specific volume than the close packed solid 
structures. The relatively dense packing of diamond 
would inhibit the phase transformations at high hy- 
drostatic pressures that are observed for heavier group 
IV elements. In addition, it suggested a revision of 
the other common notion that the diamond melting 
curve should have negative slope, something that is 
to be expected when a liquid is denser than the co- 
existing solid. Yin and Cohen also found that, at a 
pressure around 2300 GPa, diamond converts to a sim- 
ple cubic (sc) phase. This work was later extended 
pni inn 112] to consider also complex tetrahedral struc- 
tures. It was found that a distorted diamond structure 
called BC-8 was stable versus diamonds at pressures 
above 1000 GPa (see Fig. [4]). 

In 1984 Shaner and coworkers [93] shock compressed 
graphite and measured the sound velocity in the ma- 
terial at shock pressures ranging from 80 to 140 GPa, 
and corresponding shock temperatures ranging from 
1500 to 5500 K. They measured velocities close to 
those of an elastic longitudinal wave in solid dia- 
mond. These velocities are much higher than those 
of a bulk wave in a carbon melt. Since no melt 
was detected at pressures and temperatures well above 
the graphite/diamond/liquid triple point, the diamond 
melting curve should, according to these results, have 
a positive slope. In 1990 Togaya [M] reported ex- 
periments in which specimens of boron-doped semi- 
conducting diamond were melted by flash-heating at 
pressures between 6 and 18 GPa: these experiments 
provided clear indications that the melting tempera- 



ture of diamond increases with pressure. 
In the same year ab-initio molecular dynamics (MD) 
studies [IS] clearly showed that, upon melting dia- 
mond at constant density, the pressure of the system 
increases. These results imply that, at the densities 
studied, the slope of the melting curve of diamond is 
positive. The shape of the diamond melting curve has 
interesting consequences for the theory of planetary 
interiors. Given our present knowledge of the phase 
diagram of carbon and the existing estimates for the 
temperatures and pressures in the interior of the outer 
planets Neptune and Uranus, as well as in the Earth 
mantle, one might conclude that in a large fraction of 
these planetary interiors the conditions are such that 
diamond should be the stable phase of carbon [96] ; di- 
amonds could then be expected to occur wherever the 
carbon concentration is sufficiently high. In section 
[S] we show that, when only homogeneous nucleation 
is considered, our modelling predict that the driving 
force for diamond nucleation is missing in giant planet 
interiors. In 1996 Grumbach and Martin [97] made a 
systematic investigation of the solid and liquid phases 
of carbon over a wide range of pressures and temper- 
atures by using ab initio MD. These authors studied 
the melting of the simple cubic and BC-8 solid phases, 
and investigated structural changes in the liquid in the 
range 400-1000 GPa. They observed that the coordi- 
nation of the liquid changes continuously from about 
four-fold to about six-fold over this pressure range. 
In 2004, Bradley et al. [HH] reported experiments on 
laser-induced shock compression of diamond up to 3000 
GPa. Through optical reflectivity measurements, they 
found for the first time direct evidences of diamond 
melting, at an estimated pressure of P = 1000 ± 
200 GPa, and temperature T = 12000 ± 4000 K. 
Fig. [4] summarizes the information about phase dia- 
gram of carbon at the time when the present research 
was started. 

3.2 Carbon phase diagram according to 
LCBOP 

Methods We performed Monte Carlo simulations on 
the LCBOPI+ model of carbon [Ml US] to estimate the 
properties of the liquid, graphite, and diamond phases 
of carbon. Coexistence curves were determined by lo- 
cating points in the P — T diagram with equal chemical 
potential for the two phase involved. For this purpose, 
we first determined the chemical potential for the liq- 
uid, graphite, and diamond at an initial state point 
(P = 10 GPa, T = 4000 K). This state point is near the 
estimated triple point :99j. Subsequently, the liquid- 
graphite, liquid-diamond, and graphite-diamond coex- 
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Figure 4: In 2004, a schematic representation of the 
phase diagram of carbon at high pressures would have 
looked more or less as indicated in this figure. Full 
curves correspond to phase boundaries for which ther- 
modynamic data are available. More recent develop- 
ments are discussed in the next section. 



istence pressures at T = 4000 K were located. In turn, 
these coexistence points served as the starting point 
for the determination of the graphite melting line, the 
diamond melting line, and the graphite-diamond co- 
existence curve, obtained by integrating the Clausius- 
Clapeyron equation (a procedure also known as Gibbs- 
Duhem integration) |100| : jp — ^^ where Av is the 
difference in specific volume, and Ah the difference in 
molar enthalpy between the two phases. 

We proceeded by first determining the Helmholtz 
free energy at a given density and temperature by ther- 
modynamic integration and subsequently calculating 
the chemical potential using the procedure described 
in Ref. [lOlj . Coexistence at a given T is found at the 
P where the different /x cross. 

For all phases, the Helmholtz free energy F* of the 
initial state point {P = 10 GPa, T = 4000 K) was de- 
termined by transforming the system into a reference 
system of known free energy F . The transforma- 
tion was imposed by changing the interaction poten- 
tial: Ux = {1- A)t7* + XU"=^. Here, C/* and C/^^ 
denote the potential energy function of the LCBOPI+ 
and the reference system, respectively. The transfor- 
mation is controlled by varying the parameter A con- 
tinuously from to 1. The free-energy change upon 
the transformation was determined by thermodynamic 
integration: 
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(6) 



The symbol (...)a denotes the ensemble average with 
the potential Ux. 

As reference system for the liquid we chose the well- 
characterized Lennard- Jones (12-6) system, whilst the 
reference system for both diamond and graphite was 
the Einstein crystal. General guidelines for these kind 
of calculations are given in [1021 llOlj , while a full de- 
scription of the strategy adopted for the present sys- 
tems is given in [5Fj . The ensemble averages needed for 
the thermodynamic integration were determined from 
Monte Carlo (MC) simulations of a 216-particle system 
in a periodically replicated simulation box. For simu- 
lations of the graphite phase, the atoms were placed 
in a periodic rectangular box with an initial edge-size 
ratio of about 1 : 1.5 : 1.7. For the liquid phase and 
diamond a periodic cubic box was used. 

From the Helmholtz free energy to the chem- 
ical potential The chemical potential /i along the 
4000 K isotherm was obtained by integrating from the 
initial state point a fit, P{p) = a+bp+cp^ , through sim- 
ulated (P, T) state points along the 4000 K isotherm. 
Here, p is the number density, and a, b, and c are fit pa- 
rameters. This yields for the chemical potential [lOlj: 



/3m(p) 
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H-^ + b^^^+b + ^i^P-p'') 



(7) 

Here, p* denotes the number density at the initial state 
point, A'^ the number of particles, and /3 — 1/k^T, with 
Ub the Boltzmann constant. 

Calculated coexistence curves The equilibrium 
densities p^/(10'^ kg/m' ) at the initial state point 
(P = 10 GPa, T = 4000 K) were 3.425 for diamond, 
2.597 for graphite, and 2.421 for the liquid. Three con- 
figurations at the equilibrium volume were then chosen 
as starting state point for the free energy calculation. 



[GPa] b [GPa nm^] c [GPa nm^] 



Liquid 

Diamond 

Graphite 



89.972 
74.809 
108.29 



-1.9654 
-3.6307 
-2.2707 



0.011 092 
0.019 102 
0.011 925 



Table 1: Parameters for the polynomial fitting of the 
4000 K isotherms of the three phases, according to: 
P{p) = a + bp + cp'^. 

The integrals related to the reference system trans- 
formation (Eq. H]) were evaluated using a 10-point 
Gauss-Legendre integration scheme. Fig. [5] shows the 
integrand {U'''^ - c/LCBOPI+^^ versus A. The smooth 
behaviour of the curves indicates that there are no spu- 
rious phase transitions upon the transformation to the 
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Figure 5: 

^LCBOPI+ 



Plots of the quantity /3/A {W'^ - 
I \ (see Eq. [S]) as a function of the coupling 



parameter A for the liquid, graphite, and diamond. On 
the left side of the horizontal axis (A = 0) is the pure 
LCBOPI+, on the right side (A = 1) is the reference 
system, i.e. the Lennard-Jones liquid for the liquid 
phase and two Einstein crystals (with different cou- 
pling constant) for graphite and diamond. The simu- 
lated ten points per phase are marked by their error 
bars, that are almost reduced to a single dash at this 
scale. 
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Figure 6: Equations of state at 4000K for the liq- 
uid, graphite, and diamond. The curves are quadratic 
polynomial fits to the data. The circles indicate the 
points, at 10 GPa, for which the Helmholtz free en- 
ergy was determined using Eq. [Sj The solid arrows 
connect coexisting (stable) points, i.e. liquid/graphite 
and graphite/diamond. The dashed arrow indicates 
the liquid/diamond coexisting point, with graphite 
metastable. 



reference system (the absence of such transitions is a 
necessary condition for using this method). At the ini- 
tial state point (P = 10 GPa, T = 4000 K), the cal- 
culated free energies(/?.F*/A) where -24.824 ± 0.006, 
-24.583 ± 0.002, and -25.137 ± 0.002, for graphite, 
diamond, and the liquid, respectively. 

Fig. [S] shows the calculated state points along the 
4000 K isotherms for the three phases, along with 
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Figure 7: Phase diagram of carbon up to 60 GPa. The 
solid right triangle, square, and diamond are the three 
coexistence points found by equating the chemical po- 
tentials at 4000 K (see text). The open right trian- 
gles, squares, and diamonds are the calculated coex- 
istence points, propagated via Gibbs-Duhem integra- 
tion. The solid circle with error bars indicates the ex- 
perimental estimate for the liquid/graphite/diamond 
triple point [721 [Ml EH]- The dashed curve is the 
experimental graphite melting curve from Ref. [73j . 
The up triangles are graphite melting state points from 
Ref. [Tlj . The crosses represent experimental estimates 
for graphite/diamond coexistence from Ref. (S^. The 
asterisk represent the theoretical graphite/diamond co- 
existence at zero kelvin, as reported in Ref. [75] . 
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Figure 8: Phase diagram of carbon at between 
and 400 GPa. The thick solid curves are our cal- 
culated phase boundaries. The dashed curve is the 
metastable prolongation of the graphite melting curve, 
from Gibbs-Duhem integration; the curve stops just 
before the simulated graphite became unstable, dis- 
playing large density fluctuations. The dashed-dotted 
curve departing from the experimental estimate for the 
triple point (sohd circle with error bar [731 HO^j 1^ ) 
is the diamond melting curve calculated in Ref. |104j 
with the Brennerl potential. The solid circle is the fi- 
nal point of the shock wave experiment of Ref. [93] at 
which diamond is not yet melted. The crosses mark the 
(metastable) liquid state points with an equal fraction 
of three and four-fold coordinated atoms. The circles 
represent state points at which the liquid freezes. 



the fitted curves. The fit parameters are listed in 
Table [TJ Employing subsequently Eq. [7| we ob- 
tained the calculated chemical potential /i along the 
4000 K isotherm for the liquid, graphite, and dia- 
mond phase. The intersections of the chemical po- 
tential curves yield the graphite/liquid coexistence at 
6.72 ± 0.60 GPa (^icL = -24.21 ± 0.10 k^T), and 
the graphite/diamond coexistence at 15.05 ± 0.30 GPa 
(/XGD = -23.01 ± 0.03 A;bT). The third intersection 
locates a diamond/liquid coexistence at = 12.75 ± 
0.20 GPa (/XDL = -23.24 ± 0.03 k^T). Even though 
both diamond and the liquid are there metastable, the 
Clausius-Clapeyron integration of the diamond melt- 
ing curve can be started at the metastable coexistence 
point at 4000K. Starting from the three coexistence 
points at 4000 K, the coexistence curves were traced 
by integrating the Clausius-Clapeyron equation using 
the trapezoidal-rule predictor-corrector scheme jlOOj . 
The new value of the coexisting P at a given T was 
taken when two iterations differed less than 0.01 GPa, 
this being the size of the single uncertainty in the calcu- 
lation of dP/dT. Typically this required 2-3 iterations. 
The calculated phase diagram in the P — T plane 
is shown in Fig. [7| for the low pressure region, and in 
Fig. [Sj for the pressures up to 400GPa. Tab.[l|Hsts the 



densities of selected points on the coexistence curves. 
The three coexistence curves meet in a triple point at 
16.4 ± 0.7 GPa and 4250 ± 10 K. 

The graphite/diamond coexistence curve agrees well 
with the experimental data. In the region near the liq- 
uid/graphite/diamond triple point that has not been 
directly probed in experiments, the graphite/diamond 
coexistence curve bends to the right, departing from 
the commonly assumed straight line. Analysis of our 
data shows this is mainly due to the rapid reduction 
with increasing pressure of the interplanar distance in 
graphite at those premelting temperature. This causes 
an enhanced increase of the density in graphite, yield- 
ing a decrease of dT/dP. 

Table[2]shows the melting enthalpy A/i™ for graphite 
and diamond. These are calculated as the difference in 
enthalpy between the solid and the liquid at coexis- 
tence. Our calculated melting enthalpies of graphite 
are significantly lower than the values of ~ 110 kJ/mol 
that were reported in recent shock-heating melting ex- 
periments [721 m] • Nonetheless our values retain the 
feature of being rather constant along the graphite 
melting curve. To our knowledge, no experimental data 
have been reported for the melting enthalpies of dia- 
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Graphite melting curve 
P [GPa] T [K] pG 


PL 


A/i,„ 


2.00 
6.70 


3800 
4000 


2.134 
2.354 


1.759 
2.098 


68.8 
66.3 


16.4 


4250 


2.623 


2.414 


64.7 


Diamond melting curve 
P [GPa] T [K] pD 


PL 


Ahm 


16.4 


4250 


3.427 


2.414 


95.9 


25.5 


4750 


3.470 


2.607 


111.5 


43.9 


5500 


3.558 


2.870 


130.8 


59.4 


6000 


3.629 


3.043 


143.9 


99.4 


7000 


3.783 


3.264 


160.5 


148.1 
263.2 


8000 
10000 


3.960 
4.286 


3.485 
3.868 


164.7 
195.3 


330.5 


11000 


4.230 


4.045 


208.1 


408.1 


12000 


4.593 


4.236 


221.7 



Table 2: Pressure (P), temperature (T), solid and 
liquid densities (p), and melting enthalpy (A/i™) 
along the melting curves. Densities are expressed in 
10'^ kg/m^, enthalpies are in [kJ/mol]. 



mond. Note, that they monotonically increase with 
temperature. 

The calculated graphite melting temperature is 
monotonically increasing with pressure and is confined 
to small temperature range around 4000 K. In con- 
trast to data inferred from experiments it shows no 
maximum and is at a somewhat lower temperature. In 
agreement with the experiments the coexistence tem- 
perature is only slowly varying with pressure. Inspec- 
tion reveals that this behavior is due to a limited vari- 
ability of the melting enthalpy, and a similar bulk mod- 
ulus for liquid and graphite yielding a volume change 
upon melting that is almost constant along the melting 
curve. 

The sign of the slope of the diamond melting curve is 
consistent with the available experimental data [93l |73] 
(see Fig. |S]). When compared to the diamond melting 
curve of the Brennerl model [104j . the LCBOPI+ di- 
amond melting curve has a steeper slope yielding sig- 
nificantly higher temperatures for the diamond melt- 
ing curve. Recently, the melting curve of diamond in 
a range up to 2000 GPa has been studied by ab initio 
MD simulations using density functional theory. Wang 
et al. |105j determined the relative stability of the dia- 
mond and liquid phase by evaluating the free energy of 
both phases. Correa et al. fl06j determined the melting 
temperature using a "two phase" simulation method, 
where the system initially consists of a liquid and a di- 
amond structure that are in contact. Subsequently the 



melting temperature is estimated by locating the tem- 
perature at which the system spontaneously evolves 
towards a liquid or a crystalline structure. In both ab 
initio MD studies it was found that the diamond melt- 
ing curve shows a maximum; around 450 GPa |106j 
or 630 GPa |105j |l|. Subsequent laser-shock experi- 
ments [107] provided data consistent with this obser- 
vation, indicating a negative melting slope most proba- 
bly in the region of 300-500 GPa. When comparing the 
LCBOPI+ diamond melting curve, that monotonically 
increases with pressure, to the ab initio MD results 
of Refs. [1051 1106] we see a significant deviation from 
200 GPa onwards. This might be attributed to an in- 
correct description of the liquid structure at high com- 
pression. Indeed, LCBOPI+ has not been vahdated 
against high density structures with coordination be- 
yond four. These are typical configuration that might 
become more dominant in the pressure region beyond 
200 GPa. 



4 Existence of a liquid— liquid 
phase transition? 

4.1 History of the LLPT near the 
graphite melting hne 

Analysis of experimental data The possibility of a 
LLPT in liquid carbon was first investigated by Kor- 
sunskaya et al. [36_ , who analysed data on the graphite 
melting curve proposed by Bundy [65 , (those data 
showed a maximum melting temperature at 6.5 GPa). 
By fitting the data from Bundy into the original two 
levels model of Kittel [34] and postulating the existence 
of two liquids, Korsunskaya et al. found the critical 
temperature Tc of the LLPT. The model is fitted with 
three points on the graphite melting curve, with the 
respective derivatives, and with the heat of melting at 
one selected pressure. The fitting procedure gives an 
estimate for the critical pressure of ~ 6.5 GPa and for 
the critical temperature of the searched transition at 
3770 K, i.e. below the melting temperature. The fitted 
value for the entropy of freezing is the same for the two 
liquids, thus implying a vertical slope (dT/dP) of the 
coexistence curve (in the metastable liquid region just 
below the critical temperature) 

On the basis of these results, the authors were able to 
calculate also the diamond melting curve: they pre- 
dicted it to have a negative slope, in accordance with 



^The difference between these two values gives a hint on the 
uncertainties related to the two different methods used for calcu- 
lating coexistence, given that the DF-MD set-up is quite similar 
in the two works 
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the commonly accepted interpretation of the experi- 
ments of Alder and Christian [7^ . Note that the slope 
of the graphite melting curve, and the slope of the dia- 
mond/graphite coexistence, as extracted from Bundy's 
data [HQI [65] , together with the densities of the phases 
obtained by fitting to the two levels model implied (via 
Clausius-Clapeyron equation) a negative slope of the 
diamond melting curve. Different values of the slopes 
of the graphite boundary curves, and of the densities 
of the phases can yield rather different slope of the di- 
amond melting curve. 

Consistent with the slope of the fitted graphite melt- 
ing curve, the low density liquid is predicted less dense 
than the coexisting graphite, and the high density liq- 
uid more dense than the coexisting graphite. The na- 
ture of the two liquids was predicted as follows: at low 
pressure graphite melts into a liquid of neutral par- 
ticles that interact predominantly through dispersion 
(London) forces. Upon increasing pressure the liquid 
would transform into a metallic close packed liquid. No 
assumptions were made on the local structure. 
A semi-empirical equation of state The modern 
discussion on the LLPT for carbon, starts with the 
elaboration of a semi-empirical equation of state for 
carbon, valid also at high P and T, by van Thiel and 
Ree [1081 199] . This equation of state was constructed 
on the basis of experimental data and electronic struc- 
ture calculations. It postulated the existence, in the 
graphite melt, of a mixture of an sp2 and an sp^ liq- 
uid. By assuming the model of pseudo-binary mixture 
for the description of the mixing of the two liquids 38] , 
Van Thiel and Ree showed that fitting their empir- 
ical equation of state to the graphite melting points 
of Bundy [65 , they predict a graphite melting curve 
that shows a maximum with a discontinuous change of 
the slope, so that a first order LLPT arises. On the 
other hand, if they fit their model to the data from 
Ref. [inS], the predicted T^ of the LLPT drops below 
the melting curve and the transition between the two 
liquids becomes continuous in the stable liquid region. 
As pointed out by Ponyatovsky [110^ the expression 
for the mixing energy of the two liquid as proposed 
by van Thiel and Ree in [1081 [^ involves two am- 
biguities. Firstly, extrapolating the coexistence curve 
between the two liquids at atmospheric pressure, the 
coexistence temperature would be T '--^ 3700 K: this 
would imply that the sp^ liquid (and the glass) would 
be more stable than the sp2 at ambient pressure up to 
very high temperatures, which is in disagreement with 
the experimental data. Furthermore, the mixing en- 
ergy is proposed to have a linear dependence on T, so 
that, when T — > 0, also the mixing energy would tend 
to zero, i.e. at zero temperature the regular solution 



would become an ideal solution. This would be rather 
unusual. 

Experimental evidence from the graphite melt- 
ing curve Using flash-heating experiments Togaya [74] 
determined the melting line of graphite and found 
a maximum in the melting curve at Pmax — 5.6 
GPa. This author fitted the six experimental points 
with two straight lines: with positive slope at pres- 
sures lower than Pmax, with negative slope at pres- 
sures higher than Pmax- The apparent discontinu- 
ity at the maximum would imply the presence of a 
triple point graphite/ low-density-liquid (LDL) / high- 
density- liquid (HDL), as a starting point of a LLPT 
coexistence curve. 

Prediction of a short-range bond-order poten- 
tial In Ref. [103j Glosli and Rcc reported a complete 
study of a LLPT simulated with the 'Brenncrl' bond- 
order potential [52] in its version with torsional inter- 
actions [lllj . The authors simulated in the canonical 
(NVT) ensemble several samples at increasing densi- 
ties at eight different temperatures. By measuring the 
pressure, they show the familiar van der Waals loop 
betraying mechanical instability over a finite density 
range. Using the Maxwell equal-area construction, the 
authors calculated the LLPT coexistence curve, ending 
in a critical point at T = 8802 K and P = 10.56 GPa. 
The lowest temperature coexistence point was calcu- 
lated at T = 5500 K and P = 2.696 GPa. The 
LDL/HDL coexistence curve should meet the graphite 
melting curve at its maximum, but unfortunately the 
Brennerl potential does not contain non bonded inter- 
actions, thus it can describe neither bulk graphite nor 
its melting curve. To overcome this deficiency, the au- 
thors devised an ingenious perturbation method. As- 
suming constant slope of the negative sloped branch 
of the graphite melting curve (the authors adopted 
the graphite melting curve measured by Togaya 74J) 
and fixing the graphite/diamond/HDL triple point at a 
value taken from the experimental literature, they give 
an estimate of the graphite/LDL/HDL triple point, at 
T = 5133 K and P = 1.88 GPa. The LDL was found to 
be mainly two-fold (sp) coordinated with a polymeric- 
like structure, while the HDL was found to be a net- 
work forming, mainly four-fold, (sps) liquid. Following 
the predictions of this bond-order potential, the sp2 
coordinated atoms would be completely avoided in the 
liquid. The authors identified the reason in the pres- 
ence of torsional interactions. In fact, the increase in 
density demands an increase in structures with higher 
coordination than the sp, which is entropically favored 
at low densities. The single bonds of the sps structures 
can freely rotate around the bond axis, while bonds 
between sp2 sites are constrained in a (almost) planar 
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geometry by the torsional interactions: this imphes a 
low entropy for a liquid dominated by sp2 sites. This 
low entropy would eventually destabilize the sp2 sites 
towards the sp^. To prove this conjecture, the authors 
calculated two relevant isotherms in the original ver- 
sion of the potential, without torsional interactions, 
finding no sign of a LLPT. Since some torsional in- 
teractions are definitely needed to mimic the double 
bond reluctancy to twist, the authors concluded that 
the LLPT predicted by the Brenner bond-order poten- 
tial with torsion is more realistic than its absence when 
torsional interactions are switched off. 
Tight binding calculations [112] showed no evidence of 
van der Waals loops at some of the temperatures ana- 
lyzed in Ref. [103j . As Glosli and Ree note, the tight 
binding model used in jll2] is strictly two-center, thus 
the torsional interactions cannot be described. 
An ab initio study of the LLPT In Ref. gS], Wu 
et al. reported a series of NVT-CPMD simulations 
at 6000 K from density 1.27 to 3.02 10^ kg/m^, in a 
range where the Brennerl potential showed the first or- 
der LLPT at the same T. No sign of a van der Waals 
loop was found: in contrast to the Brennerl results of 
the previous paragraph, two approaching series start- 
ing from the lowest and the highest density, were found 
to meet smoothly at intermediate densities. Looking 
for the reasons of the failure of the Brennerl potential, 
the authors calculated, with the same functional used 
in the CPMD simulations, the torsional energy of two 
model molecules. One, (CH3)2C=C(CH3)2, was cho- 
sen so that the bond between the two central atoms 
represents a double bond in a carbon network: two 
sp2 sites are bonded each to two sp^ sites; the periph- 
eral hydrogens are needed to saturated the sp3 atoms 
and are intended to have no effect on the central bond. 
The second molecule, (CH2)2C-C(CH2)2 is a portion 
of a completely sp2 coordinated network: in the bond- 
order language, the central bond is conjugated. The 
two molecules were geometrically optimized in their 
planar configurations and then twisted around the cen- 
tral bond axis in steps of 7r/12. In each configuration 
the electronic wave function was optimized, without 
further relaxations, to give the total energy, that was 
compared to the planar configuration total energy. The 
difference is the torsional energy. The DF calculations 
found a surprising picture: while the double bond tor- 
sional energy was only slightly overestimated by the 
Brennerl potential at intermediate angles, the DF tor- 
sional energy for the conjugated bond showed a com- 
pletely different scenario compared to the classical pre- 
diction. It shows a maximum at 7r/4, while the planar 
and orthogonal configuration have basically the same 
energy. For the Brennerl potential, the torsional en- 



ergy in this conjugated configuration is monotonically 
increasing with the torsion angle, just as for the dou- 
ble bond configuration. On average, considering that 
the conjugated configuration would be characteristic 
of a mainly sp2 coordinated liquid, the torsional inter- 
actions are enormously overestimated by the classical 
potential. As a further illustration, the authors tried 
to lower torsional energy of the conjugated bond in 
the classical potential, by tuning the proper parame- 
ter, and found a much less pronounced LLPT. Note 
that the functional form of the torsional interactions 
for the Brennerl potential cannot reproduce the DF 
data mentioned here. Wu et al. concluded that "[the] 
Brenner potential significantly overestimates the tor- 
sional barrier of a chemical bond between two- and 
three-center-coordinatcd carbon atoms due to the in- 
ability of the potential to describe lone pair electrons" ; 
and: "[the] Brenner potential parameters derived from 
isolated hydrocarbon molecules and used in the liter- 
ature to simulate various carbon systems may not be 
adequate to use for condensed phases, especially so in 
the presence of lone pair electrons" . In the next section 
we show that the conclusion of Wu et al. is not neces- 
sarily true for all BOPs; indeed, LCBOPI+, the carbon 
bond-order potential proposed by Los and Fasolino (see 
section [2]), includes a definition of the torsional inter- 
action which is able to reproduce relevant features of 
liquid carbon, as they are described by DF-MD. 

4.2 Ruling out the LLPT in the stable 
liquid region via LCBOP 

We have already indicated that the change of the struc- 
ture of the liquid along the graphite and diamond melt- 
ing curve is related to the slope of the melting curve. 
More importantly, it plays also a crucial role in the nu- 
cleation of diamond in liquid carbon. The latter will 
be further discussed in the next section. 

The calculated melting curves of the LCBOPI"'" 
model for carbon up to 400 GPa provide strong evi- 
dence that there is no LLPT in the stable liquid phase. 
One indication is the smoothness of the slopes of the 
melting curves. A further argument lies in the struc- 
ture of the liquid near freezing. Below we discuss this 
in more detail. 

The calculated phase diagram (Figs. [7] and El) does 
not show the sharp maximum in the graphite melting 
line that was inferred from the calculated first-order 
LLPT for the Brennerl bond-order potential |103j . As 
we mentioned in the previous section, subsequent DF- 
MD simulations of liquid carbon [45l |46j indicate that 
the Brennerl LLPT is spurious: it originates from an 
inadequate description of the torsional contribution to 
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the interactions. We have extended the calculation of 
the graphite melting curve of LCBOPI+ towards higher 
pressures into the region where both graphite and the 
liquid are metastable with respect to diamond. It is 
plotted as a dashed curve in Fig. [8] that shows the same 
trend as at lower pressures. Hence, the calculated slope 
of the graphite melting curve is incompatible with the 
existence of a LLPT in this region of the carbon phase 
diagram. 

In order to further analyze the nature of the liquid, 
we determined several structural properties of the liq- 
uid near the melting curve where we also explored the 
diamond melting curve. Fig. [HI shows the coordina- 
tion fraction in the liquid along the coexistence curves 
up to 400 GPa, as function of temperature, pressure 
and density, with a linear scale in density. The dashed 
curve is the calculated graphite/diamond/liquid triple 
point. Along the graphite melting curve, the three- 
fold and two-fold coordination fractions remain rather 
constant, with the four-fold coordination slightly in- 
creasing to account for the increase in density. Along 
the diamond melting curve the three-fold coordinated 
atoms are gradually replaced by four-fold coordinated 
atoms. However, only at (3.9 10^ kg/m^, 300 GPa, 
and 10500 K) the liquid has an equal fraction of three- 
fold and four-fold coordinated atoms. The change of 
dominant coordination is rather smooth. Moreover, we 
have verified that it is fully reversible showing no sign 
of hysteresis in the region around the swapping of dom- 
inant coordination. We note, that these results contra- 
dict the generally assumed picture (see e.g. Ref. [5U] ) 
that diamond melts into a four-fold coordinated liq- 
uid. Our calculations suggest that up to ^300 GPa 
the three-fold coordination dominates. 

The interrelation between three and four-fold sites, 
was further investigated calculating the partial radial 
distribution functions {gij{r)) of the liquid at 300 GPa, 
and 10500 K. Partial radial distribution functions are 
defined as the probability of finding a j-fold site at a 
distance r from a i-fold site; the total radial distribu- 
tion function g is recovered hy: g = J2i 9ii + 2 J^i^j 9ij ■ 
We show the results in Fig. [TOl we focus on the 
three predominant curves, describing the pair corre- 
lations between three- fold atoms (533), between four- 
fold atoms (344), and the cross pair correlation be- 
tween three- and four- fold sites (534). Disregarding the 
rather pronounced minimum in correspondence of the 
dip around 2 A of the 533 and the 334, the similar- 
ity of three curves at all distances r is striking. The 
two sites are almost undistinguishable: in case of a 
tendency towards a phase transition, one would expect 
some segregation of the two structures. In contrast, 
looking at distances within the first neighbours shell, a 
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Figure 9: Coordination fraction of the liquid along 
the melting curve of carbon. The melting curves are 
unimodal, thus fixing p providing one-to-one relation 
among p, P and T. The scale for p is chosen to be lin- 
ear. The dashed curve is the liquid/graphite/diamond 
triple point. On the left hand side of the triple point, 
the liquid coexists with graphite, while on the right 
hand side it coexists with diamond. 




Figure 10: Partial distribution functions gij of the liq- 
uid at the calculated coexistence with diamond, at 
10500 K and - 300 GPa, when three- and four-fold 
atoms are equally present. The left panel is for the 
diagonal contributions (i.e. for i = j), while the right 
panel is for the cross correlations (i.e. for i 7^ j). 



three-fold site seems to bond indifferently to a three- or 
a four-fold site, and viceversa. Furthermore, the par- 
tial structures up to the third, quite pronounced, peak 
at ^ 4.5 A, are almost the same for these three partial 
radial distribution functions. 

We determined the properties of the metastable liq- 
uid in the stable diamond region. Fig.[8]shows the liq- 
uid P — T state points (crosses) that exhibit an equal 
number of three and four-fold coordinated atoms. It 
ranges from the high-pressure high-temperature region 
where the liquid is thermodynamically stable down into 
the diamond region, where the liquid is metastable for 
the LCBOP. The circles indicate state points in which 
the LCBOP liquid freezes in the simulation. Enclosed 
by the two set of points lies what we baptized diamond- 
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like liquid. This is a mainly four-fold coordinated liq- 
uid with a rather pronounced diamond-like structure in 
the first coordination shell and was discussed in [i^ . 
This suggests that a (meta) stable liquid with a domi- 
nantly four-fold coordination may only exist for pres- 
sures beyond ~ 100 GPa and could imply that the 
freezing of liquid into a diamond structure might be 
severely hindered for a large range of pressures beyond 
the graphite/diamond/liquid triple point. In Ref. [H] 
it is also pointed out that at 6000 K the equation of 
state shows a change of slope around the transition to 
the four-fold liquid. At even lower temperature this 
feature becomes more and more evident, but for tem- 
peratures lower than ~4500 K the liquid freezes into 
a mainly four-fold coordinated amorphous structure. 
This observation is consistent with quenching MD sim- 
ulations |1131I114] to obtain the tetrahedral amorphous 
carbon. In those simulations a mainly three-fold liquid 
freezes into an almost completely four-fold amorphous. 
Recent fully ab-initio study of the diamond melting 
line [1051 I106J predicted a maximum at pressures be- 
yond the maximum pressure (400 GPa) we explored 
with our potential. The maximum implies a liquid 
denser than diamond at pressure higher than the pres- 
sure at which the maximum appears. The authors of 
both works analyzed the structure of the liquid around 
this maximum, finding no sign of abrupt change in den- 
sity and/or coordination. This points towards exclud- 
ing a LLPT between a four-fold and a higher-fold coor- 
dinate liquid; rather, a smooth transformation towards 
a denser liquid is always observed. 

5 Diamond nucleation 

Our knowledge of the phase diagram of "LCBOPI"'" 
carbon" allows us to identify the regions of the phase 
diagram where diamond nucleation may occur. We 
studied the homogeneous nucleation of diamond from 
bulk liquid, by computing the steady-state nucleation 
rate and analyzing the pathways to diamond formation. 
On the basis of our calculations, we speculate that the 
mechanism for nucleation control is relevant for crys- 
tallization in many network-forming liquids, and also 
estimate the conditions under which homogeneous dia- 
mond nucleation is likely in carbon-rich stars and plan- 
ets such as Uranus and Neptune. 

Steady-state nucleation rate Most liquids can 
be cooled considerably below their equilibrium freez- 
ing point before crystals start to form spontaneously 
in the bulk. This is caused by the fact that micro- 
scopic crystallites are thermodynamically less stable 
than the bulk solid. Spontaneous crystal growth can 
only proceed when, due to some rare fluctuation, one 



or more micro-crystallites exceed a critical size (the 
"critical cluster"): this phenomenon is called homoge- 
neous nucleation. An estimate of the free-energy bar- 
rier the system has to cross in order to form critical 
clusters and of the rate at which those clusters form in 
a bulk super-cooled liquid, can be obtained from Clas- 
sical Nucleation Theory (CNT) (nil- CNT assumes 
that AG(n), the Gibbs free-energy difference between 
the metastable liquid containing an n-particle crystal 
cluster and the pure liquid, is given by 



AG{n) = S{n)j-n\Anl 



(8) 



where S{n) is the area of the interface between an n- 
particle crystallite and the metastable liquid, 7 is the 
liquid-solid surface free-energy per unit area, and Afi 
the difference in chemical potential between the solid 
and the super-cooled liquid. The surface area S{n) is 
proportional to c{n/ps)'^^^, where the factor c depends 
on the shape and the geometry of the cluster (e.g. c = 
167r/3 for a spherical cluster). 

The top of the free-energy barrier AG* to grow the 
crystalline critical cluster is then given by 



AG* 



r 
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(9) 



where ps is the number density of the stable phase and 
c indicates the geometrical properties of the growing 
cluster. From our simulations, we can only determine 
the product 07: it is this quantity and the degree of 
super-saturation (A/x), that are needed to compute the 
top of the free-energy barrier, and hence the nucleation 
rate. 

CNT relates R, the steady-state nucleation rate, i.e. 
the number of crystal clusters that form per second 
per cubic meter, to AG*, the height of the free-energy 
barrier that has to be crossed to nucleate the critical 
crystal n*: 



R 



CNT 



-/3AG(n*) 



(10) 



where AG* is the top of the free-energy barrier and k 
is the kinetic prefactor. The kinetic prefactor term is 
defined as 

K^ PLk+^n'Z (11) 

where pL is the liquid number density, /c+^n* the attach- 
ment rate of single particles to a spherical crystalline 
cluster k+^n' = (24D(7i*)2/3) /A^, with D/A^ propor- 
tional to the jump frequency (A being the atomic jump 
distance) and Z — y^\Ap\/{6TrkBTn*) the so-called 
Zeldovitch factor. As the nucleation rate depends ex- 
ponentially on AG*, a doubling of 7 may change the 
nucleation rate by many orders of magnitude. 
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Because of the extreme conditions under which ho- 
mogeneous diamond nucleation takes place, there have 
been no quantitative experimental studies to determine 
its rate. Moreover, there exist no numerical estimates 
of A/x and 7 for diamond in super-cooled liquid car- 
bon. Hence, it was thus far impossible to make even 
an order-of-magnitude estimate of the rate of diamond 
nucleation. 

Results We simulate a 2744 particles bulk liquid 
carbon using periodic boundary conditions, with a cu- 
bic box whose edge is 18 A. We make it metastable 
by undercooling it at constant pressure at two differ- 
ent state points, A{P = 85 GPa, T = 5000 K} and 
B{P = 30 GPa, T = 3750 K}. At both state points, 
the liquid is super-cooled by {Tm — T)/Tm ~ 25% be- 
low the diamond melting curve, Tm being the melting 
temperature T^ = 6600 K and T^ = 5000 K, respec- 
tively. 

We evaluate A/i by thermodynamic integration at 
constant pressure from the melting point (/3m = 
I/ZcbTm) 



A(AAi) 



13. 



{[hsW)-hLm)pdp (12) 
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where f3i — l/fcsTi, {i = A, B}, and h is the enthalpy 
per particle of the solid and liquid phase, respectively. 
We then find: \AnA/kBT\ = 0.60 and lA/xs/fesTl = 
0.77. 

In recent years several authors have been developing 
methods for studying homogeneous nucleation from the 
bulk and detecting solid particles within the metastable 
liquid [nnmnilllllini]- our study on diamond nu- 
cleation is based on these works, but it requires var- 
ious adaptations due to the specificity of the carbon 
covalent bond [55] . We have already shown that liquid 
carbon is rather structured below its freezing curve. 
This leads to the need of building a " strict" definition 
of "crystallinity" of a particle, in order to avoid an 
overestimation of the number of solid particles in the 
system. 

In order to compute the nucleation free energy, we 
use the biggest crystal cluster n as a local order param- 
eter to quantify the transformation from the liquid to 
the solid. To identify solid-like particles, we analyze the 
local environment of a particle using a criterion based 
on a spherical-harmonics expansion of the local bond 
order. In practice, the present bond-order parameter is 
based on rotational invariants constructed out of rank 
three spherical harmonics (Ysm)- This choice allows 
us to identify the tetragonal symmetry of the diamond 
structure, as already described in Ref. [5511119) . and it 
is also perfectly suited to find particles in a graphite- 
like environment. Our choice of odd-order of spherical 



harmonics is due to the fact that both diamond and 
graphite lattices have odd symmetry upon inversion of 
coordinates. 

In order to define the local order parameter, we start 
with computing 
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93,; 






(13) 



where the sum extends over all neighbors of particle i 
and over all values of m. Zi is the fractional number of 
neighbours and S'^°'^^(rij) is a smooth cut-off function, 
introduced in the context of LCBOP/^ [55] (see also 
Section II. A in [iU]). 

By properly normalizing Eq. I13[ we get 
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being q^ ^ the complex conjugate of 93,™ ■ 

Next we define the dot product between the normal- 
ized function q^ ^ of particle i and the same function 
computed for each of its first neighbors, d3{i,j), and 
sum them up over all the m values: 
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d,{i,3) = E '?3,™(0 • 9;%(j)^""""(r-.,)- (15) 



d'i{i^ j) is a real number defined between -1 and 1: it as- 
sumes the value of -1 when computed for both graphite 
and diamond ideal structures. 

Two neighboring particles i and j are considered 
to be connected whenever d3{i,j) < dc — —0.87. 
This value satisfactorily splits the distributions of solid 
particles belonging to a thermalized lattice and liq- 
uid particles as found in a liquid. The histograms 
that led us to this choice are thoroughly discussed 
in Refs. (551 1119] . By counting the total number of 
connections (jicon) and plotting the probability distri- 
bution of ricon, we define a threshold for the number 
of connections needed to neatly distinguish between 
a liquid- like and a solid- like environment: we assume 
that whenever ricon > i^con = 3 8' particle is solid-like. 
At this stage, we do not specify any nature of the par- 
ticle's crystallinity, whether diamond-like or graphite- 
like. By means of a cluster algorithm we then define 
all the solid-like AND connected particles as belong- 
ing to the same crystal cluster. After computing the 
size of each cluster, we use the size of the biggest clus- 
ter as the order parameter which describes the phase 
transition [120j . 

Once properly identified the biggest crystalline clus- 
ter in the system, we use the umbrella sampling tech- 
nique [121] to measure the free-energy barrier AG* to 
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form a critical cluster at state point A. In order to 
better equilibrate the growing clusters, we implement 
a "parallel tempering" algorithm similar to the one de- 
scribed in Ref. |122j . We obtain that, at state point 
A, AG\ is around 25 ksT for a critical cluster size of 
UA — 110. By fitting the initial slope of Eq. [5]to a poly- 
nomial function assuming a spherical growing cluster, 
while imposing the value of the correspondent super- 
saturation {PA/i =0.60), the inter-facial free energy is 
7A = 0.21kBT/k^ ~ 1.86 J m-2. The same value of 
7^ is obtained from the top of the free-energy barrier 
assuming a spherical cluster shape (Eq. [5]) [ps = 0.191 
A""^). We underline the fact that at the chosen ther- 
modynamic conditions, there are no finite size effects, 
caused by spurious interaction of the critical cluster 
with its own periodically repeated image. 

By knowing the inter-facial free energy in A, and 
assuming the validity of CNT, we estimate the crys- 
tal nucleation rate by means of Eq. [TOl where we use 
Eq. [TT] to compute the kinetic pre-factor (the atomic 
jump distance A being of the order of the diamond bond 
distance, 1.54A): i?^^^'^ - ©(lO^O) s'^m-^. 

We also use Forward-Flux Sampling (FFS), a rela- 
tively recent rare events technique useful to compute 
the nucleation rate and to study the pathways to nu- 
cleation |123lll24|lll9j . and we measure the crystal nu- 
cleation rate at state point A. FFS yields an estimate 
for the nucleation rate that is three orders of magni- 
tude higher than the one estimated by means of Eq. 1101 
Whilst such a discrepancy seems large, it need not be 
significant because nucleation rates are extremely sen- 
sitive to small errors in the calculation of the nucleation 
barrier. Two possible reasons for this discrepancy are: 

1) if we consider that the standard deviation corre- 
sponding to 7 is around 10% of its measured value, we 
conclude that the nucleation rate is 0(10^°^^) s~^m~'^; 

2) another source of error can be the poor statistics 
when computing the nucleation rate from molten car- 
bon by means of FFS. This is due to the time con- 
suming calculations of the interaction potential: in our 
study we are in fact forced to base our results on 0(10) 
independent nucleation events. Fig. Illl shows a typical 
critical cluster at state point A obtained in the FFS 
simulations: it contains around 110 particles, and it is 
surrounded by mainly 4-fold coordinated liquid parti- 
cles. The picture shows two different views of the same 
cluster: it appears evident that all particles within the 
bulk are diamond-like, whereas the particles belonging 
to the outer surface are less connected but still mainly 
3-4 fold coordinated. 

We then attempt to compute the nucleation rate at 
state point B by means of FFS and, even in rather 
long simulations, we cannot observe the formation of 








^f>«^■^■:i^/^' 



Figure 11: Two diflterent views of the biggest cluster at 
state point A containing around 110 particles, surrounded 
by mainly 4-fold coordinated liquid particles. 



any crystal cluster containing more than 75 particles. 
Hence, these calculations suggest that the nucleation 
rate at state point B measured by means of FFS is 
around zero. Figurefl^shows a 75 particles cluster con- 




Figure 12: Typical snapshot of a crystalline cluster of ^ 75 
particles obtained at state point B, surrounded by mainly 
3-fold coordinated liquid particles. 



taining 3-fold coordinated surface particles surround- 
ing the 4-fold coordinated bulk particles, while embed- 
ded in a 2-3 fold coordinated liquid. 

As we are unable to grow critical nuclei with FFS, 
we assume that a system of 2744 particles is too small 
to accommodate a spherical critical cluster. According 
to Classical Nucleation Theory (CNT) [TT5] . the crys- 
tal nucleation rate depends exponentially on the height 
of the free-energy barrier (see Eq. [T(I| . The latter is a 
function of the inter-facial free energy (7) cube and 
inversely proportional to the super-saturation (A/x) 
square. Since the super-saturation is quite similar in 
both state points, the failure of the system to nucleate 
suggests that the inter-facial free energy should play a 
major role. In order to estimate the free-energy barrier 
in state point B, as we know the solid number density 
{ps = 0.177 A~'^) and the chemical potential differ- 
ence between the liquid and the solid {(3Ap.B =0.77), 
we only need to calculate the inter-facial free energy 7. 
Thus, in what follows, we focus on methods to estimate 
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7 at state point B. 

As a spherical critical cluster does not fit in our sim- 
ulation box, we prepare a rod-like crystal in a system 
with a slab geometry: this is a flattened box contain- 
ing around 4000 particles, with lateral dimensions that 
are some four times larger than its height. The crystal 
rod is oriented perpendicular to the plane of the slab, 
it spans the height of the simulation box and is con- 
tinued periodically. The cross section of this crystal 
rod is initially lozenge shaped, such that its [lll]-faces 
are in contact with the liquid. The [lll]-planes are 
the most stable ones for the diamond lattice. In fact, 
macroscopic natural diamonds have often an octahe- 
dral shape, with eight [lll]-exposed surfaces. Indeed, 
we find stable [111] surfaces in all but the smallest stud- 
ied diamond clusters. 

At state point A clusters grow by the addition of 
particles to the surface made of mainly [lll]-planes. 
Note that when graphite and diamond structures com- 
pete at state point B (as shown in Fig. [T2|) . the [0001]- 
graphite sheets transform into [111] diamond planes. 
Fig. [T3] represents the top view at state point _B of a 




Figure 13: Top view of a rectangular parallelepiped formed 
by 4 [lll]-faces and 2 [l-10]-lozenge bases with an acute 
angle 0=70.52 degrees thermalized at state point B. 



rod- like crystal, formed by 4 [lll]-faces and 2 bases as 
[l-10]-lozenge with the acute angles of 0=70.52 degrees. 
We then rewrite Eq. \E\ for a rectangular parallelepiped 
having 4 faces and 2 lozenge-shaped basis 
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where h is the slab's height. We then use Umbrella 
Sampling to compute the initial slope of the free- 
energy barrier. As ft. = 10 A, we obtain from fit- 
ting Eq. [in] that the inter-facial free energy for the 
lozenge-shaped cluster is 7^ = 0.91kBT/k'^ ~ 4.70 



J m~^. At the same time, computing the inter- facial 
free energy of the rod-like crystal at state point A gives 
7^ = 0.37fcBT/A^ ~ 2.55 J m~^, considering the same 
slab's height and the same angle 9. 

Now that we have estimates for the inter-facial free 
energies of the lozenge-shaped clusters at both state 
points A and B, we can estimate the ratio between 
them and find that c^b/ctia = 1b Ma — 2.5. As we 
compare clusters having the same shape, this ratio is 
presumably not very sensitive to the precise (and, a 
priori unknown) shape of the cluster shape. As the sur- 
face free energy at state point B is appreciably higher 
than at state point A, the early stages of crystal forma- 
tion at point B are strongly suppressed by the inter- 
facial free-energy term. Since we know 7^ (referred 
to a hypothetical spherical cluster) and the ratio be- 
tween the two 7's, we can infer that the "effective" 
7s = 0.68fcsr/A2 ~ 3.50 J m-^. 

This value of the surface free energy is so large that 
we would indeed have needed a much larger system in 
order to accommodate the critical cluster at the state 
B thermodynamic conditions. /,From Eq. 1161 we cal- 
culate the critical cluster size for the lozenge-shaped 
parallelepiped ri^Y^ and use it to estimate the size of a 
critical spherical cluster rigp in B\ 
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At state point B we find B ri^Y^ ~330 particles. Ex- 
pressing ngp as a function of the lozenge-shaped par- 
allelepiped one, we get 
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where ps is the solid number density ps = 0.17 A ^, 
[A^B/fesTl = 0.77, and h the height of the slab (10 A). 
Thus, njD ~ 1700 particles at state point B. To guar- 
antee that the critical cluster does not interact with its 
own periodic images, its radius should always be less 
than 25% of the box diameter L. A spherical cluster 
with a radius of 0.25 L occupies ~ 7 % of the volume 
of the box and, as the solid is denser than the liquid, 
it contains about 10 % of the total number of parti- 
cles {N w 17000). Such a large system size is beyond 
our present computational capacity. In contrast, in the 
slab geometry we find that the free energy of a lozenge- 
shaped crystal goes through a maximum at a size of ~ 
330 particles, which is much less than the system size 
(4000 particles). 

As A/i and ps are known, we can now use CNT to es- 
timate AG* in state point B. It turns out that, mainly 
because jb is 2.5 times larger than ja, the nucleation 
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barrier in B is more than ten times higher than in point 
A, and the nucleation rate is Rb ^ 10^^" s~^m~^. 

To understand the microscopic origin for the large 
difference in nucleation rates in state points A and B, 
it is useful to compare the local structure of the liquid 
phase in both state points. As discussed in section Wl\ 
above (see also [13 [US]), liquid carbon is mainly 4- fold 
coordinated at state point A (20% 3-fold and 80% 4- 
fold ), whereas at the lower temperatures and pressures 
of point B, the coordination in the liquid resembles 
that of graphite and is mainly 3-fold coordinated ( 5% 
2-fold, 85% 3-fold and 10% 4-fold ). 

We can analyze the structure of the crystalline clus- 
ters that form in the supersaturated liquid carbon and 
distinguish graphite-like from diamond-like particles. 
In an a posteriori analysis, we use a different order pa- 
rameter function of the order two spherical harmonics, 
and particularly sensitive to the graphite planar ge- 
ometry. q2m{i) is the linear combination of spherical 
harmonics computed for each particle i 
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where the sum extends over all neighbors of particle 
i. We then sum over all the m values and calculate 
the modulus, \q2\- The \q2\ probability distribution for 
both A and B is represented in Fig.[T31 Figure [TJ] de- 
picts the features of both the smallest (~ 20 in both 
state points A and B) and the biggest clusters {^ 250 
in A and ~ 75 in B). We also distinguish among: 
liquid- like particles (circles) , particles belonging to the 
surface of the largest cluster (squares), particles inside 
the bulk cluster (diamonds) and particles belonging to 
the first liquid layer surrounding the largest cluster (tri- 
angles). According to our definition, particles belong- 
ing to the surface of the cluster are those connected to 
solid-like particles, but not solid-like themselves. Con- 
cerning particles belonging to the first liquid layer sur- 
rounding the cluster, they usually display the same be- 
haviour as the ones belonging to the cluster surface, 
which is not surprising in view of the uncertainty in dis- 
tinguishing a surface-particle from a first-liquid-layer 
particle. To neatly distinguish between diamond-like 
or graphite- like environment, we use as a reference the 
1 52 1 probability distribution for both bulk diamond (D) 
and graphite (G) (inset of Fig. [T4|). 

At state point A, it is clear that bulk particles be- 
longing to small clusters (bottom-left side) and big 
clusters (bottom-right side) are mainly diamond-like, 
as well as particles belonging to the surface of the 
clusters. In contrast, at state point B bulk particles 
belonging to small clusters (top-left side) show both 
graphite-like and diamond-like finger-prints. By vi- 
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Figure 14: The top-left represents clusters of ^ 20 and 
the top-right clusters of ~ 75 at state point B, whereas 
the bottom-left clusters of ^ 20 and the bottom-right clus- 
ters of ~ 250 particles at state point A. The used code 
is: circles=liquid particles, squares=particles belonging to 
the biggest cluster's surface, diamonds=particles within the 
bulk cluster, triangles=particles belonging to the first liq- 
uid layer surrounding the biggest cluster. The inset shows 
the \q2\ probability distribution for an equilibrated bulk di- 
amond (D) (left-hand side) and graphite (G) (right-hand 
side) . 



sual inspection, we note that when clusters grow larger 
(around 75 particles), particles at the surface tend to 
be mainly 3-fold coordinated, whereas bulk particles 
stay 4-fold coordinated, as shown in the top-right side 
of Fig. [T31 The destabilizing effect of the graphitic liq- 
uid on the diamond clusters is most pronounced for 
small clusters (large surface-to- volume ratio): clusters 
containing less than 25 particles tend to be graphitic in 
structure, clusters containing up to 60 particles show a 
mixed graphite-diamond structure (see Fig. [T^ . It ap- 
pears that the unusual surface structure of the diamond 
cluster is an indication of the poor match between a di- 
amond lattice and a 3-fold coordinated liquid. 

5.1 Consequences for other network 
forming liquids, carbon-rich stars, 
Uranus and Neptune 

As discussed in the introduction, there are many 
network- forming liquids that, upon changing pressure 
and temperature, undergo profound structural changes 
or even LLPT [551 [23 I126J . Interestingly, our simula- 
tions show that the ease of homogeneous crystal nucle- 
ation at constant super-saturation from one-and-the- 
same meta-stable liquid can be tuned by changing its 
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pressure, and thereby its local structure. 

Pressures and temperatures that we investigate for 
the diamond nucleation are in practice impossible to 
reach in experiments. However, such conditions are 
likely to be found in several extraterrestrial "laborato- 
ries" . Homogeneous nucleation of diamond may have 
taken place in the atmosphere of carbon-rich binary 
stellar systems comprising the so-called carbon stars 
and white dwarfs [1271 1128] . Closer to home, it has 
been suggested that diamonds could also have formed 
in the carbon-rich middle layer of Uranus and Nep- 
tune [Ml 11291 I130J where, due to the high pressure 
and temperature, the relatively abundant CH4 would 
decompose into its atomic components. In fact, ex- 
periments on methane laser-heated in diamond anvil 
cells [131) found evidence for diamond production. Ab 
initio simulations [132j also found that hot, compressed 
methane will dissociate to form diamond. Yet, there is 
a large discrepancy between the estimates of the pres- 
sures (and thus depth in the planet interior) at which 
the diamond formation would take place. The laser- 
heating experiments jl31j suggested diamond forma- 
tion at pressure as low as 10-20 GPa (at 2000-3000 K), 
whereas the ab initio simulations |132j found dissocia- 
tion of methane, but synthesis of short alkane-chains at 
~ 100 GPa and diamond at pressures not lower than 
300 GPa (note that simulations were carried out at 
4000-5000 K). 

The present work allows us to make a rough estimate of 
the conditions that are necessary to yield appreciable 
diamond nucleation on astronomical timescales. 

In this context, it is crucial to note that neither car- 
bon stars nor carbon-rich planets consist of pure car- 
bon. In practice, the carbon concentration may be as 
high as ~50% in carbon- rich stars |1271ll28j . but much 
less (1-2% [HHllMlIinailSn]) in Uranus and Neptune. 
To give a reference point, it is useful to estimate an up- 
per bound to the diamond nucleation rate by consider- 
ing the rate at which diamonds would form in a hypo- 
thetical environment of pure, metastable liquid carbon. 
To this end we use our numerical data on the chem- 
ical potential of liquid carbon and diamond and our 
numerical estimate of the diamond-liquid surface free 
energy, to estimate the nucleation barrier of diamond 
as a function of temperature and pressure. We then 
use CNT to estimate the rate of diamond nucleation. 

To do so, we need to extend the estimate of the nu- 
cleation rate from the triple point pressure (around 16 
GPa) up to 100 GPa, and from the melting tempera- 
tures ( T^ = 6600 K and T,^ = 5000 K, respectively) 
to 35 % under-cooling (at which diffusion in our sam- 
ple becomes negligible on the - far from astronomical 
- time-scales of our simulations). To make such an 



extrapolation, we make use of Eqn.'s [TOl [H and [TT] 
The state-point dependent quantities are the solid and 
liquid number densities pL and ps, the self-diffusion co- 
efficient Z?, the surface free energy 7, the difference in 
chemical potential between the liquid and the solid Ap, 
and the critical cluster size n*. We estimate them in 
the following way: the densities are directly measured 
by Monte Carlo simulations of the solid and the liquid; 
the self diffusion coefficient is extrapolated assuming 
an Arrhenius behaviour of the metastable liquid (see 
Appendix B); the chemical potential difference is in- 
terpolated via Eq. [T^] between 30 and 85 GPa. Ap 
then also follows by linear interpolation. Concerning 
the surface free-energy, we assume that 7(P, T) linearly 
depends on C4, the equilibrium concentration of 4- fold 
coordinated atoms at the selected state point. This 
quantity is easily measured in the Monte Carlo simula- 
tions. The nucleation barrier height is given by Eq. 
where the geometrical factor c is the same for all cluster 
sizes. It is obvious that we have to make rather drastic 
assumptions in order to estimate the nucleation rate 
in the experimentally relevant regime. We believe that 
our assumptions are reasonable, but one should not 
expect the resulting numbers to provide more than a 
rough indication. 
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Figure 15: The figure shows part of the carbon phase di- 
agram from Ref. [134] and the iso-nucleation rate zones. 
The solid red curves represent the coexistence curves from 
Ref. [m]. Pa = 85 GPa, Ta = 5000 K and Pb = 30 GPa, 
Tb ~ 3750 K. Along the green dashed curve the ratio of 
3- fold and 4- fold coordination in the liquid is 1:1. The 
numbers on the right indicate the base 10 logarithm (or 
the order of magnitude) of the crystal nucleation rate from 
molten carbon (in m~^s~^). The continuous black curve 
is the boundary of the region above which nucleation rate 
becomes negligible (< 10~'*''m^'^s~^). 
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Figure dni shows that there is a region of some lOOOK 
below the freezing curve (continuous red curve) where 
diamond nucleation is less than 10^"*° m^^s^^ (above 
the continuous black line) . If the rate is lower than this 
number, not a single diamond could have nucleated in 
a Uranus-sized body during the life of the universe. 
As can be seen from the figure, our simulations for 
state point B are outside the regime where observable 
nucleation would be expected. Note that this latter 
conclusion is not based on any extrapolation. 

As mentioned above, carbon stars and planets do not 
consist of pure carbon. Hence, we have to consider the 
effect of dilution on the crystallization process. To do 
so, we make a very "conservative" assumption, namely 
that nucleation takes place from an ideal mixture of C, 
N,0 and H [135j . If this were not the case, then either 
demixing would occur, in which case we are back to the 
previous case, or the chemical potential of carbon in the 
liquid is lower than that in pure carbon, which would 
imply that the thermodynamic driving force for dia- 
mond crystallization is less than in pure liquid carbon. 
In Fig. 1161 we show how dilution affects the regime 
where diamond nucleation is possible. To simplify this 
figure, we do not vary pressure and temperature inde- 
pendently but assume that they follow the adiabatic 
relation that is supposed to hold along the isentrope of 
Uranus |136j and we use the ideal-mixture expression 
for the chemical potential (3Afi = /3A^o + /3^"-([C])j 
where /3A/zo is the chemical potential difference be- 
tween the solid and the liquid for he pure substance 
(C) and [C] is the concentration of carbon in the fluid 
mixture. 

Not surprisingly. Fig. [11] shows that dilution of the 
liquid decreases the driving force for crystallization. In 
fact, no stable diamond phase is expected for carbon 
concentrations below 8%. Moreover, there is a wide 
range of conditions where diamonds could form in prin- 
ciple, but never will in practice. Assuming that, for a 
given pressure, the width of this region is the same as in 
the pure C case (almost certainly a serious underesti- 
mate), we arrive at the estimate in Fig.[Tl]of the region 
where nucleation is negligible (i.e. less than one dia- 
mond per planet per life-of-the-universe). From this 
figure, we see that quite high carbon concentrations 
(over 15%) are needed to get homogeneous diamond 
nucleation. Such conditions do exist in white dwarfs, 
but certainly not in Uranus or Neptune. 



Appendix A: LCBOPII 

In this appendix, we describe the main features of the 
latest addition to the LCBOP family. This potential 
has been used in Ref. jiUl [5H H^T] . However, the sim- 
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Figure 16: Diamond nucleation boundary as a function of 
carbon concentration: in the plot, the rate is zero (no ther- 
modynamic driving force to nucleation) in the top region 
(liquid), it is negligible (< 10~*°m~'^s~^) in the middle 
region and non-negligible (> 10~*°7Ti~^s~^) in the bottom- 
right region. We call the nucleation rate negligible if it cor- 
responds to less than one cluster per Uranus-sized planet 
over a period of 10^" years. The left hand y-axis rep- 
resents the temperature; the right-hand y-axis indicates 
the corresponding pressure for a Uranus-like isentrope (see 
Ref. [T29ll96l [T33llT30| l. 



Illations discussed in the present paper are based on 
LCBOPI+. 

Middle range interactions. Although LCBOFI"*" 
gave an improved description of most liquid phase 
properties, like coordination distributions as a function 
of density, as compared to the bond-order potentials 
without LR interactions (Brenner, CBOP P^B] ). the ra- 
dial distribution function showed a too marked mini- 
mum after the first shell of neighbors, as compared to 
ab-initio calculations (see Fig. [^3]). This deficiency was 
attributed to the relatively short cut-off of 2.2 A for 
the SR interactions, giving rise to a spurious barrier for 
bond formation around 2.1 A. Therefore, for LCBOPII, 
the total binding energy expression was extended with 
so-called MR interactions as: 



E,^-y { S^V^ + (1 - 5^[) Vl; + ^^5""U"'^ 

9 / -/ \ '-J '■J V ''J / ''J nyrnr '-J '-J 

(20) 

The first two terms on the right-hand side represent 
the SR and LR interactions respectively, where S'f'" 
smoothly switches between both interactions within 
the interval 1.7 A< r^ < 2.2 A, with S'^''(1.7) = 1 
and S"'''(2.2) = 0. The last term represents the MR in- 
teractions, where Uj™'' is a purely attractive potential 
and Z™'' is a sort of MR coordination number defined 
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Figure 17: Schematic representation of the domains of 
the three types of interactions. The numbers and 1 
below the hnes indicate the values of the switch func- 
tions as well as the intervals where these switch func- 
tions are applied. As an example, the SR interactions 
is smoothly switched off between rij — 1.7 A and 
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tion S^U"- 



to account for many body effects. The switch func- 
'S'""'(?'y), going from to 1 between rij 

1.7 A and r^ 2.2 A, smoothly excludes the MR 
interactions for distances smaller than 1.7 A. For clar- 
ity, the ranges of the various interactions in Eq. [^D] 
are schematically represented in Fig. [T71 The MR 
interaction was fitted to ab-initio calculations of sin- 
gle, double and triple bond dissociation curves. For 
the single bond, the tail of the interaction vanishes be- 
yond 4 A. Vj^^ is the product of a simple polynomial 
VJPf^ with a smooth cut-off at 4 A and an environment 
dependent switch function Sg^^{{9ijk}), depending on 
the angles between r.^ — Yj — r^ and r^^ = r^ — r^ 
(Vfc ^ j), where atom fc is a SR neighbour of atom z, 
i.e. \Yik\ < 2.2. Thus, while the MR interactions give 
an extension of the covalent interactions beyond the 
SR cut-off distance 2.2 A in situations where this is 
appropriate, its environment dependence relies only on 
the SR nearest neighbours (within 2.2 A), a quite con- 
venient property for the sake of efficiency. The switch 
S'fP'^ acts in such a way that V"j^^ is only non-zero when 
the angles 9ijk are relatively large. This is illustrated 
schematically in Fig. 1181 In particular, the definition 
of S™[, makes V^Y^ vanish for any pair ij in all bulk 
crystal structures. So the addition of the MR interac- 
tions does not require reparametrization of the SR and 
LR potential terms. 

The reactivity of atoms depends on whether these 
atoms are well surrounded by neighbours or not. Typ- 
ically, an atom with a dangling bond wants to make 
another bond. To include this effect, the MR poten- 



tial is made dependent on the so-called dangling bond 
number Nf\ i.e. V^[^ = V^Jj{Nf). For an atom with 
Nf^ ~ 1, the MR interaction is stronger than for an 
atom with Nf'' = 0. 



Extended coordination dependence of angu- 
lar function. For LCBOPII the correction of the an- 
gle dependent part of the bond order for configurations 
involving low coordinations and small angles has been 
further extended, involving a gradual coordination de- 
pendence of the angular term over a wide range of co- 
ordinations. 

Anti-bonding. Another new feature of LCBOPII is 
the addition of an anti-bonding correction to the bond 
order. An example of a situation where one electron 
remains unpaired in a non-bonding state is depicted 
in Fig. fTglb). Clearly, this situation is unfavourable 
as compared to the situation in Fig. flW c). This ef- 
fect cannot be captured in the conjugation term and 
has therefore been included as a separate, anti-bonding 
term. 

Torsion. As it has been clearly demonstrated in 
Ref. |45j, the torsion interaction for a bond ij be- 
tween an atom i and an atom j is strongly dependent 
on conjugation, i.e. on the coordinations of the neigh- 
bours k{^ j) of i and l{y^ i) of j. This was already 
partly included in LCBOPI+. However, for LCBOPII, 
the conjugation dependence of the torsion interactions 
was fully extended and fitted to ab-initio calculations 
of the torsion barrier for all the possible conjugation 
situations. 

In addition to that, LCBOPII includes a redefini- 
tion of the torsion angle, in order to avoid the 'spuri- 
ous' torsion that occurs using the traditional definition. 




Figure 18: Example of a typical situations where 
MR interactions are active (upper graph) and where 
they are switched off by the switch function Sg (lower 
graph), since Sg vanishes for small bond angles Ojik 
and Oiji. 
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Figure 19: Illustration of configurations where the 
bond energy is not monotonously dependent on A^™""* , 
due to the fact that bonding is less effective when the 
electron supply from atom i, Nf^, is not equal to that 
from atom j, N^^. A correct description of each of the 
three cases was achieved by adding a negative anti- 
bonding term to the bond order which depends on 
AiVj'^- — Nfj — Nj-. The resulting binding energy is 
given below each configuration (see Fig. [^0)1 . 



Traditionally, the torsion angle LUijki is defined as 



COs{uJijkl) 



t'ijk 



^ijk I l^ijl I 



tiji _ (r,j X r^fc) ■ (r^j x r^i) 
|ry X rikWvij X Tjil 



(22) 

which, assuming torsion to be non-vanishing only be- 
tween sp^ bonded atoms i and j, gives rise to four 
torsion contributions to the bond order. With the def- 
inition of Eq. [221 both situations depicted in Fig. [2Tb 
and c give rise to a non-zero torsion angle. However, in 
both situations there is actually no torsional distortion 
but only a bending distortion, which is already taken 
into account by the angular term in the bond order. 
Thus, one would like to have cJijki — for the cases in 
Fig. 121b and c, in disagreement with the most right- 
hand side expression in Eq. [551 Another problem of 
expression 1221 is that it has a singularity for configura- 
tions where r^j is parallel to r,;fc (or ru ) . For the Hquid 
phase at high temperature such situation are easily ac- 
cessible. 

For LCBOPII, the problem of 'spurious' torsion has 
been tackled by a redefinition of the vectors tijk in Eq. 
[551 reading: 



^ijk 



X (fifcj -fifes) + 



V3 



+ — (i"'U • {i^k^ - f.fej) (f,j X (f,fc, + f.fcj) 

and likewise for tju. Inserting these vectors into Eq. 
[551 leading to a different right-hand side, reproduces 
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Figure 20: Torsional barriers according to LCBOPII 
and DF calculations for the six possible values of N^^^-' 
for an ij bond between sp^ atoms i and j. Sym- 
bols represent the DF results, curves the fits obtained 
by the LCBOPII. Top panel: torsional barriers for 
the extreme values of N^°^\ related to the conju- 
gated (A^'^°"^=0, squares and dashed curve) and dou- 
ble bonds (iV'^°"-' = l, circles and solid curve). Bot- 
tom panel: intermediate values of N'^°"'^: 1/8 (stars 
and dotted curve), 1/4 (down triangles and dashed- 
dotted curve), 1/2 (up triangles and solid curve) and 
5/8 (diamonds and dashed curve). Note the complex 
behaviour of the curves for the values 1/2 and 5/8, 
where the barrier at tt/2 is higher for iV^°"-'=l/2 than 
for N'^°"^=5/8. The top panel alone applies also to 
LCBOPI+ (see Appendix A in :40j for details.) 




(23) Figure 21: Illustration of the occurrence of spurious 
torsion when using the definition of Eq. [221 (most right- 
hand side) for the torsion angles. 
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the same Uijki as the traditional definition for any tor- 
sional distortion without bending and yields tOijki = 
for both situations depicted in Fig. HU as it should be. 
In addition, it gives a good interpolation for any other 
configuration, and it has no singularities. Note that 
for the two distortions depicted in Fig. [^ the second 
term in Eq. 1241 vanishes, and the vectors tijk and tju 
are parallel, implying indeed ujijki — 0. 



configurations with integer configurations, the weight 
factors Wc for the configuration c (=1,..4) being defined 
in terms of the switch functions Sz{rik)- For instance, 
for the situation in Fig. 1221 the conjugation term is 
given by: 




Figure 22: For LCBOPII, a fractional coordination sit- 
uation is treated as a weighted superpositions of integer 
coordination situations, as illustrated in this picture. 
Dashed lines indicate partial bonds, solid lines are for 
full bonds. 

Interpolation for fractional coordinations. The 

conjugation term F^°^^\Nij,Nji,N^°^-') for a bond ij 
depends on the reduced coordinations TVy and Nji of 
the atoms i and j, and on the conjugation number 
iV™"-'. The reduced coordination Nij is defined as: 






(24) 



and likewise for Nji, where Sz(rik) is a switch func- 
tion for the coordination, smoothly going from 1 to 
for nk going from 1.7 A and 2.2 A. F[^"^ is fit- 
ted to integer coordination configurations with only 
full neighbours, i.e. with Sz{rik) = Sz{rji) = 1 ^k,l. 
This poses the problem of how to determine i^™""' for 
configurations with fractional bonds, i.e. configura- 
tions with Sz{rik) < 1 for one or more neighbours k. 
David Brenner, the inventor of the conjugation term, 
proposed to use a 3D spline [52]. However, since the 
values on the integer argument nodes are rather scat- 
tered, a spline unavoidably introduces unphysical os- 
cillations. For LCBOPII, we found an alternative solu- 
tion to this problem which is schematically illustrated 
in Fig. 1221 In this approach, the conjugation term for 
configurations with fractional coordination is defined 
as a weighted superposition of conjugation terms for 
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where N^°"'' are the conjugation numbers for the four 
configurations in Fig. [22l 

Results with LCBOPII. LCBOPII proved to be 
more accurate than its predecessors in describing de- 
fects and surfaces of the solid phases [101 ■ ^^ the liquid 
phases the improvement of LCBOPII is immediately 
evident when looking at radial distribution functions at 
different densities (Fig. [23|l . The main discrepancy be- 
tween LCBOPI"*" and the reference data from DF-MD 
calculations [461 I54j ) was found at the first minimum, 
at around 2 A. LCBOPI"'" predicted a much deeper 
minimum than DF-MD. This discrepancy is completed 
eliminated by LCBOPII. We also know that the melt- 
ing line of diamond predicted by LCBOPII is about 
500 K lower than for LCBOPI+ at - 60 GPa (Sij, 
more consistently with ab-initio predictions of the dia- 
mond melting line [10511106] . In [H] we also thoroughly 
analyzed the properties of the liquid. Interestingly, 
by extrapolating the equations of state of the liquid 
at temperatures at which our relatively small samples 
actually froze, we found that a critical point for the 
graphite-like into diamond-like transition is present at 
1230 K. The precise value might be inaccurate, since it 
is found far outside the sampled region, still the shapes 
of the higher temperature equations of state point to- 
wards the existence of such critical isotherm. As is the 
case for the much speculated water LLPT [^ [5U] an 
unreachable critical point might still be responsible of 
some peculiar behaviour of the system at higher tem- 
peratures, such as the enormous change in nucleation 
rate with pressure. 

Appendix B: Self-diffusion coefficient 

When computing the kinetic pre-factor to get the nu- 
cleation rate, we have to consider the fact that for 
our model potential, LCBOPI+, only a Monte Carlo 
code is available. In order to evaluate the self-diffusion 
coefficient needed to compute the CNT kinetic pre- 
factor, we infer the scaling factor between the Monte 
Carlo "time-step" and the MD time-step [139] by prop- 
agating a 128 carbon atoms system via Car-Parrinello 
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Figure 23: Comparison of the radial distribution func- 
tions at 6000 K and four selected densities between 
LCBOPII, LCBOPI+, and the reference data taken 
from our own DF-MD simulations. 
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Figure 24: Arrhenius plot for the super-cooled liquid car- 
bon. The activation energy is Ea = 7eV from Ref. [14l| . 



Fedosayev [142] reported a measurement of the molten 
carbon viscosity: 77 = 5 xlO^^ poise at T = ISGOiiT. We 
estimate the self-diffusion coefficient at the same tem- 
perature by means of the Stokes-Einstein relation |143j 



D = 



ksT 



(25) 



where a — 1.54^ and ks is the Boltzmann's constant: 
D{m60K) = 3.3 X 10-" cmVs. Since molten carbon 
is an Arrhenius- like fluid |144| . 



D{T)^ Doexp "B-^ 



(26) 



we obtain Dq: Dq — 470cm^/s, and then extrapolate 
the diffusion coefRcient for different temperatures, as 
shown in Fig. [M] and Table [31 



T[K] 


D[AVs] 


T[K] 


D[AV.] 


1750 


2.1 X 10"^ 


3000 


6.3 X 10" 


2080 


3.5 X 10^ 


3300 


7.6 X 10^ 


2380 


5.0 X 10* 


3570 


5.0 X 10» 


2680 


2.4 X 10^ 


3760 


2.0 X 10^ 



Table 3: Self diffusion coefficient as a function of temper- 
ature. 

We then find that at state point A, Da — 3.5 x 10"^ 
cm^/s, whereas at state point B, Db = 2x 10"^ cm^/s. 

We also use a Car-Parrinello Molecular Dynam- 
ics [140] to calculate the self-diffusion coefficient by 
means of the mean square displacement: at state point 
A D = 2.3 X 10-^cm^/s, which matches surprisingly 
well with the diffusion coefficient estimated by means 
of the Arrhenius law, D = 3.5 x 10"^ cm^/s. 



Molecular Dynamics CPMD code 140 starting from 
a configuration equilibrated with LCBOPI"*" . Note the 
reasonably good agreement between the static proper- 
ties of the liquid carbon computed with LCBOPI+ and 
the same computed by means of CPMD (with the BP 
functional) [J^l [M] . Data for the high pressure state 
point come from simulations used in Ref. 46J , whereas 
data for the low pressure state point come from a new 
simulation where the time rescaling is state-point de- 
pendent, obtained with the same technical details as 
reported in [^ . 

We use the fact that molten carbon is an Arrhenius- 
like liquid: therefore, once the activation energy is 
known, we compute the viscosity as a function of tem- 
perature and by means of the Stokes-Einstein relation 
obtain the diffusion coefficient. In the nineteen-fifties, 
Kanter [141j estimated the relevant activation energy 
of liquid carbon to be Ea =683 -^. Subsequenty, 
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